Changeset 4797 for palm/trunk/SOURCE/ocean_mod.f90
- Timestamp:
- Nov 26, 2020 4:02:39 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/ocean_mod.f90 (modified) (112 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/ocean_mod.f90
r4768 r4797 1 1 !> @file ocean_mod.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 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/>. 5 ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General 6 ! Public License as published by the Free Software Foundation, either version 3 of the License, or 7 ! (at your option) any later version. 8 ! 9 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the 10 ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 ! Public License for more details. 12 ! 13 ! You should have received a copy of the GNU General Public License along with PALM. If not, see 14 ! <http://www.gnu.org/licenses/>. 16 15 ! 17 16 ! Copyright 2017-2020 Leibniz Universitaet Hannover 18 !-------------------------------------------------------------------------------- !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4768 2020-11-02 19:11:23Z suehring 27 29 ! Enable 3D data output also with 64-bit precision 28 30 ! … … 32 34 ! 4671 2020-09-09 20:27:58Z pavelkrc 33 35 ! Implementation of downward facing USM and LSM surfaces 34 ! 36 ! 35 37 ! 4535 2020-05-15 12:07:23Z raasch 36 38 ! bugfix for restart data format query 37 ! 39 ! 38 40 ! 4517 2020-05-03 14:29:30Z raasch 39 41 ! added restart with MPI-IO for reading local arrays 40 ! 42 ! 41 43 ! 4495 2020-04-13 20:11:20Z raasch 42 44 ! restart data handling with MPI-IO added 43 ! 45 ! 44 46 ! 4481 2020-03-31 18:55:54Z maronga 45 47 ! vector directives added to force vectorization on Intel19 compiler 46 ! 48 ! 47 49 ! 4346 2019-12-18 11:55:56Z motisi 48 ! Introduction of wall_flags_total_0, which currently sets bits based on static 49 ! topographyinformation used in wall_flags_static_050 ! 50 ! Introduction of wall_flags_total_0, which currently sets bits based on static topography 51 ! information used in wall_flags_static_0 52 ! 51 53 ! 4329 2019-12-10 15:46:36Z motisi 52 54 ! Renamed wall_flags_0 to wall_flags_static_0 53 ! 55 ! 54 56 ! 4272 2019-10-23 15:18:57Z schwenkel 55 ! Further modularization of boundary conditions: moved boundary conditions to 56 ! respective modules 57 ! Further modularization of boundary conditions: moved boundary conditions to respective modules 57 58 ! 58 59 ! 4196 2019-08-29 11:02:06Z gronemeier 59 60 ! Consider rotation of model domain for calculating the Stokes drift 60 ! 61 ! 61 62 ! 4182 2019-08-22 15:20:23Z scharf 62 63 ! Corrected "Former revisions" section 63 ! 64 ! 64 65 ! 4110 2019-07-22 17:05:21Z suehring 65 ! Pass integer flag array as well as boundary flags to WS scalar advection 66 ! routine 67 ! 66 ! Pass integer flag array as well as boundary flags to WS scalar advection routine 67 ! 68 68 ! 4109 2019-07-22 17:00:34Z suehring 69 69 ! implemented ocean_actions 70 ! 70 ! 71 71 ! 3767 2019-02-27 08:18:02Z raasch 72 ! unused variable for file index and tmp_2d removed from rrd-subroutine parameter 73 ! list 74 ! 72 ! unused variable for file index and tmp_2d removed from rrd-subroutine parameter list 73 ! 75 74 ! 3719 2019-02-06 13:10:18Z kanani 76 ! Changed log_point to log_point_s, otherwise this overlaps with 77 ! 'all progn.equations'cpu measurement.78 ! 75 ! Changed log_point to log_point_s, otherwise this overlaps with 'all progn.equations' 76 ! cpu measurement. 77 ! 79 78 ! 3684 2019-01-20 20:20:58Z knoop 80 79 ! nopointer option removed 80 ! 81 81 ! 3294 2018-10-01 02:37:10Z raasch 82 82 ! initial revision … … 89 89 ! Description: 90 90 ! ------------ 91 !> This module contains relevant code for PALM's ocean mode, e.g. the 92 !> prognostic equation for salinity, the equation of state for seawater,93 !> the Craik Leibovich force (Stokes force), and wavebreaking effects94 !------------------------------------------------------------------------------ !91 !> This module contains relevant code for PALM's ocean mode, e.g. the prognostic equation for 92 !> salinity, the equation of state for seawater, the Craik Leibovich force (Stokes force), and wave 93 !> breaking effects 94 !--------------------------------------------------------------------------------------------------! 95 95 MODULE ocean_mod 96 97 98 USE arrays_3d, &99 ONLY: prho, prho_1, rho_ocean, rho_1, sa, sa_init, sa_1, sa_2, sa_3,&100 sa_ p, tsa_m, flux_l_sa, flux_s_sa, diss_l_sa, diss_s_sa101 102 USE control_parameters, &103 ONLY: atmos_ocean_sign, bottom_salinityflux, 104 constant_top_salinityflux, restart_data_format_output, ocean_mode, top_salinityflux,&105 w all_salinityflux, loop_optimization, ws_scheme_sca96 97 98 USE arrays_3d, & 99 ONLY: diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, prho, prho_1, rho_1, rho_ocean, sa, & 100 sa_1, sa_2, sa_3, sa_init, sa_p, tsa_m 101 102 USE control_parameters, & 103 ONLY: atmos_ocean_sign, bottom_salinityflux, constant_top_salinityflux, loop_optimization,& 104 ocean_mode, restart_data_format_output, top_salinityflux, wall_salinityflux, & 105 ws_scheme_sca 106 106 107 107 USE kinds 108 108 109 USE pegrid, &109 USE pegrid, & 110 110 ONLY: threads_per_task 111 111 112 USE statistics, &112 USE statistics, & 113 113 ONLY: sums_wssas_ws_l 114 114 115 USE indices, &116 ONLY: advc_flags_s, n xl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0, nbgp115 USE indices, & 116 ONLY: advc_flags_s, nbgp, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0 117 117 118 118 USE restart_data_mpi_io_mod, & … … 120 120 wrd_mpi_io_global_array 121 121 122 USE surface_mod, & 123 ONLY: bc_h, surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, & 124 surf_usm_h, surf_usm_v 122 USE surface_mod, & 123 ONLY: bc_h, surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v 125 124 126 125 USE exchange_horiz_mod, & … … 132 131 CHARACTER (LEN=20) :: bc_sa_t = 'neumann' !< namelist parameter 133 132 134 INTEGER(iwp) :: ibc_sa_t !< integer flag for bc_sa_t 135 INTEGER(iwp) :: iran_ocean = -1234567 !< random number used for wave breaking effects 136 133 INTEGER(iwp) :: ibc_sa_t !< integer flag for bc_sa_t 134 INTEGER(iwp) :: iran_ocean = -1234567 !< random number used for wave breaking effects 137 135 INTEGER(iwp) :: sa_vertical_gradient_level_ind(10) = -9999 !< grid index values of sa_vertical_gradient_level(s) 138 136 139 137 LOGICAL :: salinity = .TRUE. !< switch for using salinity 140 138 LOGICAL :: stokes_force = .FALSE. !< switch to switch on the Stokes force 139 LOGICAL :: surface_cooling_switched_off = .FALSE. !< variable to check if surface heat flux has been switched off 141 140 LOGICAL :: wave_breaking = .FALSE. !< switch to switch on wave breaking effects 142 LOGICAL :: surface_cooling_switched_off = .FALSE. !< variable to check if surface heat flux has been switched off143 141 144 142 REAL(wp) :: alpha_wave_breaking = 3.0_wp !< coefficient for wave breaking generated turbulence from Noh et al. (2004), JPO … … 153 151 REAL(wp) :: u_star_wave_breaking !< to store the absolute value of friction velocity at the ocean surface 154 152 155 REAL(wp), DIMENSION(12), PARAMETER :: nom = &156 (/ 9.9984085444849347D2, 7.3471625860981584D0, &157 -5.3211231792841769D-2, 3.6492439109814549D-4, &158 2.5880571023991390D0, -6.7168282786692354D-3, &159 1.9203202055760151D-3, 1.1798263740430364D-2, &160 9.8920219266399117D-8, 4.6996642771754730D-6, &161 -2.5862187075154352D-8, -3.2921414007960662D-12 /)162 !< constants used in equation of state for seawater163 164 REAL(wp), DIMENSION(13), PARAMETER :: den = &165 (/ 1.0D0, 7.2815210113327091D-3, &166 -4.4787265461983921D-5, 3.3851002965802430D-7, &167 1.3651202389758572D-10, 1.7632126669040377D-3, &168 -8.8066583251206474D-6, -1.8832689434804897D-10, &169 5.7463776745432097D-6, 1.4716275472242334D-9, &170 6.7103246285651894D-6, -2.4461698007024582D-17, &171 -9.1534417604289062D-18 /)172 !< constants used in equation of state for seawater153 REAL(wp), DIMENSION(12), PARAMETER :: nom = & 154 (/ 9.9984085444849347D2, 7.3471625860981584D0, & 155 -5.3211231792841769D-2, 3.6492439109814549D-4, & 156 2.5880571023991390D0, -6.7168282786692354D-3, & 157 1.9203202055760151D-3, 1.1798263740430364D-2, & 158 9.8920219266399117D-8, 4.6996642771754730D-6, & 159 -2.5862187075154352D-8, -3.2921414007960662D-12 /) 160 !< constants used in equation of state for seawater 161 162 REAL(wp), DIMENSION(13), PARAMETER :: den = & 163 (/ 1.0D0, 7.2815210113327091D-3, & 164 -4.4787265461983921D-5, 3.3851002965802430D-7, & 165 1.3651202389758572D-10, 1.7632126669040377D-3, & 166 -8.8066583251206474D-6, -1.8832689434804897D-10, & 167 5.7463776745432097D-6, 1.4716275472242334D-9, & 168 6.7103246285651894D-6, -2.4461698007024582D-17, & 169 -9.1534417604289062D-18 /) 170 !< constants used in equation of state for seawater 173 171 174 172 SAVE 175 173 176 PUBLIC :: bc_sa_t, ibc_sa_t, prho_reference, sa_surface, 177 sa_vertical_gradient , sa_vertical_gradient_level,&178 sa_vertical_gradient_level_ind, stokes_force,wave_breaking174 PUBLIC :: bc_sa_t, ibc_sa_t, prho_reference, sa_surface, sa_vertical_gradient, & 175 sa_vertical_gradient_level, sa_vertical_gradient_level_ind, stokes_force, & 176 wave_breaking 179 177 180 178 … … 285 283 ! 286 284 !-- Add INTERFACES that must be available to other modules (alphabetical order) 287 PUBLIC eqn_state_seawater, ocean_actions, ocean_check_data_output, & 288 ocean_check_data_output_pr, ocean_check_parameters, & 289 ocean_data_output_2d, ocean_data_output_3d, ocean_exchange_horiz, & 290 ocean_define_netcdf_grid, ocean_header, ocean_init, & 291 ocean_init_arrays, ocean_parin, ocean_prognostic_equations, & 292 ocean_swap_timelevel, ocean_rrd_global, ocean_rrd_local, & 293 ocean_wrd_global, ocean_wrd_local, ocean_3d_data_averaging, & 294 ocean_boundary_conditions, stokes_drift_terms, wave_breaking_term 285 PUBLIC eqn_state_seawater, ocean_actions, ocean_check_data_output, ocean_check_data_output_pr, & 286 ocean_check_parameters, ocean_data_output_2d, ocean_data_output_3d, & 287 ocean_exchange_horiz, ocean_define_netcdf_grid, ocean_header, ocean_init, & 288 ocean_init_arrays, ocean_parin, ocean_prognostic_equations, ocean_swap_timelevel, & 289 ocean_rrd_global, ocean_rrd_local, ocean_wrd_global, ocean_wrd_local, & 290 ocean_3d_data_averaging, ocean_boundary_conditions, stokes_drift_terms, & 291 wave_breaking_term 295 292 296 293 297 294 CONTAINS 298 295 299 !------------------------------------------------------------------------------! 300 ! Description: 301 ! ------------ 302 !> Equation of state for seawater as a function of potential temperature, 303 !> salinity, and pressure. 296 !--------------------------------------------------------------------------------------------------! 297 ! Description: 298 ! ------------ 299 !> Equation of state for seawater as a function of potential temperature, salinity, and pressure. 304 300 !> For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech. 305 301 !> eqn_state_seawater calculates the potential density referred at hyp(0). 306 302 !> eqn_state_seawater_func calculates density. 307 303 !> TODO: so far, routine is not adjusted to use topography 308 !------------------------------------------------------------------------------ !304 !--------------------------------------------------------------------------------------------------! 309 305 SUBROUTINE eqn_state_seawater 310 306 311 USE arrays_3d, &307 USE arrays_3d, & 312 308 ONLY: hyp, prho, pt_p, rho_ocean, sa_p 313 USE indices, &309 USE indices, & 314 310 ONLY: nxl, nxr, nyn, nys, nzb, nzt 315 311 … … 355 351 sa2 = sa1 * sa1 356 352 357 pnom = nom(1) + nom(2)*pt1 + nom(3)*pt2 + &358 nom(4)*pt3 + nom(5)*sa1 + nom(6)*sa1*pt1 + &353 pnom = nom(1) + nom(2)*pt1 + nom(3)*pt2 + & 354 nom(4)*pt3 + nom(5)*sa1 + nom(6)*sa1*pt1 + & 359 355 nom(7)*sa2 360 356 361 pden = den(1) + den(2)*pt1 + den(3)*pt2 + &362 den(4)*pt3 + den(5)*pt4 + den(6)*sa1 + &363 den(7)*sa1*pt1 + den(8)*sa1*pt3 + den(9)*sa15 + &357 pden = den(1) + den(2)*pt1 + den(3)*pt2 + & 358 den(4)*pt3 + den(5)*pt4 + den(6)*sa1 + & 359 den(7)*sa1*pt1 + den(8)*sa1*pt3 + den(9)*sa15 + & 364 360 den(10)*sa15*pt2 365 361 ! … … 367 363 prho(k,j,i) = pnom / pden 368 364 369 pnom = pnom + nom(8)*p1 + nom(9)*p1*pt2 + &365 pnom = pnom + nom(8)*p1 + nom(9)*p1*pt2 + & 370 366 nom(10)*p1*sa1 + nom(11)*p2 + nom(12)*p2*pt2 371 367 372 pden = pden + den(11)*p1 + den(12)*p2*pt3 + &368 pden = pden + den(11)*p1 + den(12)*p2*pt3 + & 373 369 den(13)*p3*pt1 374 370 … … 380 376 381 377 ! 382 !-- Neumann conditions are assumed at top boundary and currently also at 383 !-- bottom boundary(see comment lines below)378 !-- Neumann conditions are assumed at top boundary and currently also at bottom boundary 379 !-- (see comment lines below) 384 380 prho(nzt+1,j,i) = prho(nzt,j,i) 385 381 rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i) … … 411 407 412 408 413 !------------------------------------------------------------------------------ !409 !--------------------------------------------------------------------------------------------------! 414 410 ! Description: 415 411 ! ------------ 416 412 !> Same as above, but for grid point i,j 417 !------------------------------------------------------------------------------ !413 !--------------------------------------------------------------------------------------------------! 418 414 SUBROUTINE eqn_state_seawater_ij( i, j ) 419 415 420 USE arrays_3d, &416 USE arrays_3d, & 421 417 ONLY: hyp, prho, pt_p, rho_ocean, sa_p 422 418 423 USE indices, &419 USE indices, & 424 420 ONLY: nzb, nzt 425 421 … … 446 442 REAL(wp) :: sa2 !< temporary scalar user for calculating density 447 443 444 448 445 DO k = nzb+1, nzt 449 446 ! … … 463 460 sa2 = sa1 * sa1 464 461 465 pnom = nom(1) + nom(2)*pt1 + nom(3)*pt2 + &462 pnom = nom(1) + nom(2)*pt1 + nom(3)*pt2 + & 466 463 nom(4)*pt3 + nom(5)*sa1 + nom(6)*sa1*pt1 + nom(7)*sa2 467 464 468 pden = den(1) + den(2)*pt1 + den(3)*pt2 + &469 den(4)*pt3 + den(5)*pt4 + den(6)*sa1 + &470 den(7)*sa1*pt1 + den(8)*sa1*pt3 + den(9)*sa15 + &465 pden = den(1) + den(2)*pt1 + den(3)*pt2 + & 466 den(4)*pt3 + den(5)*pt4 + den(6)*sa1 + & 467 den(7)*sa1*pt1 + den(8)*sa1*pt3 + den(9)*sa15 + & 471 468 den(10)*sa15*pt2 472 469 ! … … 474 471 prho(k,j,i) = pnom / pden 475 472 476 pnom = pnom + nom(8)*p1 + nom(9)*p1*pt2 +&473 pnom = pnom + nom(8)*p1 + nom(9)*p1*pt2 + & 477 474 nom(10)*p1*sa1 + nom(11)*p2 + nom(12)*p2*pt2 478 pden = pden + den(11)*p1 + den(12)*p2*pt3 + & 475 476 pden = pden + den(11)*p1 + den(12)*p2*pt3 + & 479 477 den(13)*p3*pt1 480 478 … … 510 508 511 509 512 !------------------------------------------------------------------------------ !510 !--------------------------------------------------------------------------------------------------! 513 511 ! Description: 514 512 ! ------------ 515 513 !> Equation of state for seawater as a function 516 !------------------------------------------------------------------------------ !514 !--------------------------------------------------------------------------------------------------! 517 515 REAL(wp) FUNCTION eqn_state_seawater_func( p, pt, sa ) 518 516 … … 532 530 REAL(wp) :: sa2 !< temporary scalar user for calculating density 533 531 532 534 533 ! 535 534 !-- Pressure is needed in dbar … … 549 548 550 549 551 eqn_state_seawater_func = &552 ( nom(1) + nom(2)*pt1 + nom(3)*pt2 + nom(4)*pt3 + &553 nom(5)*sa + nom(6)*sa*pt1 + nom(7)*sa2 + nom(8)*p1 + &554 nom(9)*p1*pt2 + nom(10)*p1*sa + nom(11)*p2 + nom(12)*p2*pt2 &555 ) / &556 ( den(1) + den(2)*pt1 + den(3)*pt2 + den(4)*pt3 + &557 den(5)*pt4 + den(6)*sa + den(7)*sa*pt1 + den(8)*sa*pt3 + &558 den(9)*sa15 + den(10)*sa15*pt2 + den(11)*p1 + den(12)*p2*pt3 + &559 den(13)*p3*pt1 &560 )550 eqn_state_seawater_func = & 551 ( nom(1) + nom(2)*pt1 + nom(3)*pt2 + nom(4)*pt3 + & 552 nom(5)*sa + nom(6)*sa*pt1 + nom(7)*sa2 + nom(8)*p1 + & 553 nom(9)*p1*pt2 + nom(10)*p1*sa + nom(11)*p2 + nom(12)*p2*pt2 & 554 ) / & 555 ( den(1) + den(2)*pt1 + den(3)*pt2 + den(4)*pt3 + & 556 den(5)*pt4 + den(6)*sa + den(7)*sa*pt1 + den(8)*sa*pt3 + & 557 den(9)*sa15 + den(10)*sa15*pt2 + den(11)*p1 + den(12)*p2*pt3 + & 558 den(13)*p3*pt1 & 559 ) 561 560 562 561 … … 564 563 565 564 566 !------------------------------------------------------------------------------ !565 !--------------------------------------------------------------------------------------------------! 567 566 ! Description: 568 567 ! ------------ 569 568 !> Reads the ocean parameters namelist 570 !------------------------------------------------------------------------------ !569 !--------------------------------------------------------------------------------------------------! 571 570 SUBROUTINE ocean_parin 572 571 … … 576 575 577 576 578 NAMELIST /ocean_parameters/ bc_sa_t, bottom_salinityflux, salinity, &579 sa_surface, sa_vertical_gradient, sa_vertical_gradient_level,&580 stokes_waveheight, stokes_wavelength, surface_cooling_spinup_time,&581 top_salinityflux, wall_salinityflux, wave_breaking577 NAMELIST /ocean_parameters/ bc_sa_t, bottom_salinityflux, salinity, sa_surface, & 578 sa_vertical_gradient, sa_vertical_gradient_level, & 579 stokes_waveheight, stokes_wavelength, surface_cooling_spinup_time,& 580 top_salinityflux, wall_salinityflux, wave_breaking 582 581 583 582 ! … … 607 606 END SUBROUTINE ocean_parin 608 607 609 !------------------------------------------------------------------------------ !608 !--------------------------------------------------------------------------------------------------! 610 609 ! Description: 611 610 ! ------------ 612 611 !> Check parameters routine for the ocean mode 613 !------------------------------------------------------------------------------ !612 !--------------------------------------------------------------------------------------------------! 614 613 SUBROUTINE ocean_check_parameters 615 614 616 USE control_parameters, & 617 ONLY: coupling_char, coupling_mode, initializing_actions, & 618 message_string, use_top_fluxes 619 620 USE pmc_interface, & 615 USE control_parameters, & 616 ONLY: coupling_char, coupling_mode, initializing_actions, message_string, use_top_fluxes 617 618 USE pmc_interface, & 621 619 ONLY: nested_run 622 620 … … 638 636 ! 639 637 !-- Check ocean setting 640 IF ( TRIM( coupling_mode ) == 'uncoupled' .AND. & 641 TRIM( coupling_char ) == '_O' .AND. & 638 IF ( TRIM( coupling_mode ) == 'uncoupled' .AND. TRIM( coupling_char ) == '_O' .AND. & 642 639 .NOT. ocean_mode ) THEN 643 640 644 641 ! 645 !-- Check whether an (uncoupled) atmospheric run has been declared as an 646 !-- ocean run (this settingis done via palmrun-option -y)647 message_string = 'ocean mode does not allow coupling_char = "' // 648 TRIM( coupling_char ) //'" set by palmrun-option "-y"'642 !-- Check whether an (uncoupled) atmospheric run has been declared as an ocean run (this setting 643 !-- is done via palmrun-option -y) 644 message_string = 'ocean mode does not allow coupling_char = "' // TRIM( coupling_char ) // & 645 '" set by palmrun-option "-y"' 649 646 CALL message( 'ocean_check_parameters', 'PA0317', 1, 2, 0, 6, 0 ) 650 647 … … 665 662 ibc_sa_t = 1 666 663 ELSE 667 message_string = 'unknown boundary condition: bc_sa_t = "' // & 668 TRIM( bc_sa_t ) // '"' 664 message_string = 'unknown boundary condition: bc_sa_t = "' // TRIM( bc_sa_t ) // '"' 669 665 CALL message( 'ocean_check_parameters', 'PA0068', 1, 2, 0, 6, 0 ) 670 666 ENDIF … … 673 669 674 670 IF ( .NOT. salinity ) THEN 675 IF ( ( bottom_salinityflux /= 0.0_wp .AND. & 676 bottom_salinityflux /= 9999999.9_wp ) .OR. & 677 ( top_salinityflux /= 0.0_wp .AND. & 678 top_salinityflux /= 9999999.9_wp ) ) & 671 IF ( ( bottom_salinityflux /= 0.0_wp .AND. bottom_salinityflux /= 9999999.9_wp ) .OR. & 672 ( top_salinityflux /= 0.0_wp .AND. top_salinityflux /= 9999999.9_wp ) ) & 679 673 THEN 680 message_string = 'salinityflux must not be set for ocean run ' // & 681 'without salinity' 674 message_string = 'salinityflux must not be set for ocean run ' // 'without salinity' 682 675 CALL message( 'ocean_check_parameters', 'PA0509', 1, 2, 0, 6, 0 ) 683 676 ENDIF … … 685 678 686 679 IF ( ibc_sa_t == 1 .AND. top_salinityflux == 9999999.9_wp ) THEN 687 message_string = 'boundary condition: bc_sa_t = "' // &688 TRIM( bc_sa_t ) //'" requires to set top_salinityflux'680 message_string = 'boundary condition: bc_sa_t = "' // TRIM( bc_sa_t ) // & 681 '" requires to set top_salinityflux' 689 682 CALL message( 'ocean_check_parameters', 'PA0069', 1, 2, 0, 6, 0 ) 690 683 ENDIF 691 684 692 685 ! 693 !-- A fixed salinity at the top implies Dirichlet boundary condition for 694 !-- salinity. In this case specification of a constant salinity flux is 695 !-- forbidden. 696 IF ( ibc_sa_t == 0 .AND. constant_top_salinityflux .AND. & 697 top_salinityflux /= 0.0_wp ) THEN 698 message_string = 'boundary condition: bc_sa_t = "' // & 699 TRIM( bc_sa_t ) // '" is not allowed with ' // & 700 'top_salinityflux /= 0.0' 686 !-- A fixed salinity at the top implies Dirichlet boundary condition for salinity. In this case 687 !-- specification of a constant salinity flux is forbidden. 688 IF ( ibc_sa_t == 0 .AND. constant_top_salinityflux .AND. top_salinityflux /= 0.0_wp ) THEN 689 message_string = 'boundary condition: bc_sa_t = "' // TRIM( bc_sa_t ) // & 690 '" is not allowed with top_salinityflux /= 0.0' 701 691 CALL message( 'ocean_check_parameters', 'PA0070', 1, 2, 0, 6, 0 ) 702 692 ENDIF … … 707 697 stokes_force = .TRUE. 708 698 ELSE 709 IF ( ( stokes_waveheight <= 0.0_wp .AND. stokes_wavelength > 0.0_wp ) & 710 .OR. & 711 ( stokes_waveheight > 0.0_wp .AND. stokes_wavelength <= 0.0_wp ) & 712 .OR. & 713 ( stokes_waveheight < 0.0_wp .AND. stokes_wavelength < 0.0_wp ) ) & 699 IF ( ( stokes_waveheight <= 0.0_wp .AND. stokes_wavelength > 0.0_wp ) .OR. & 700 ( stokes_waveheight > 0.0_wp .AND. stokes_wavelength <= 0.0_wp ) .OR. & 701 ( stokes_waveheight < 0.0_wp .AND. stokes_wavelength < 0.0_wp ) ) & 714 702 THEN 715 message_string = 'wrong settings for stokes_wavelength and/or ' // & 716 'stokes_waveheight' 703 message_string = 'wrong settings for stokes_wavelength and/or stokes_waveheight' 717 704 CALL message( 'ocean_check_parameters', 'PA0460', 1, 2, 0, 6, 0 ) 718 705 ENDIF … … 722 709 723 710 724 !------------------------------------------------------------------------------ !711 !--------------------------------------------------------------------------------------------------! 725 712 ! Description: 726 713 ! ------------ 727 714 !> Check data output. 728 !------------------------------------------------------------------------------ !715 !--------------------------------------------------------------------------------------------------! 729 716 SUBROUTINE ocean_check_data_output( var, unit ) 730 717 731 718 IMPLICIT NONE 732 719 … … 751 738 752 739 753 !------------------------------------------------------------------------------ !740 !--------------------------------------------------------------------------------------------------! 754 741 ! Description: 755 742 ! ------------ 756 743 !> Check data output of profiles 757 !------------------------------------------------------------------------------ !744 !--------------------------------------------------------------------------------------------------! 758 745 SUBROUTINE ocean_check_data_output_pr( variable, var_count, unit, dopr_unit ) 759 746 760 USE arrays_3d, &747 USE arrays_3d, & 761 748 ONLY: zu, zw 762 749 763 USE control_parameters, &750 USE control_parameters, & 764 751 ONLY: data_output_pr 765 752 … … 772 759 IMPLICIT NONE 773 760 761 CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit 774 762 CHARACTER (LEN=*) :: unit !< 775 763 CHARACTER (LEN=*) :: variable !< 776 CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit777 764 778 765 INTEGER(iwp) :: var_count !< 766 779 767 780 768 SELECT CASE ( TRIM( variable ) ) … … 837 825 838 826 839 !------------------------------------------------------------------------------ !827 !--------------------------------------------------------------------------------------------------! 840 828 ! Description: 841 829 ! ------------ 842 830 !> Define appropriate grid for netcdf variables. 843 831 !> It is called out from subroutine netcdf. 844 !------------------------------------------------------------------------------ !832 !--------------------------------------------------------------------------------------------------! 845 833 SUBROUTINE ocean_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) 846 834 847 835 IMPLICIT NONE 848 836 … … 854 842 LOGICAL, INTENT(OUT) :: found !< flag if output variable is found 855 843 844 856 845 found = .TRUE. 857 846 … … 860 849 SELECT CASE ( TRIM( var ) ) 861 850 862 CASE ( 'rho_sea_water', 'rho_sea_water_xy', 'rho_sea_water_xz', &863 ' rho_sea_water_yz', 'sa', 'sa_xy', 'sa_xz', 'sa_yz' )851 CASE ( 'rho_sea_water', 'rho_sea_water_xy', 'rho_sea_water_xz', 'rho_sea_water_yz', 'sa', & 852 'sa_xy', 'sa_xz', 'sa_yz' ) 864 853 grid_x = 'x' 865 854 grid_y = 'y' … … 877 866 878 867 879 !------------------------------------------------------------------------------ !868 !--------------------------------------------------------------------------------------------------! 880 869 ! Description: 881 870 ! ------------ 882 871 !> Average 3D data. 883 !------------------------------------------------------------------------------ !872 !--------------------------------------------------------------------------------------------------! 884 873 SUBROUTINE ocean_3d_data_averaging( mode, variable ) 885 886 887 USE arrays_3d, &874 875 876 USE arrays_3d, & 888 877 ONLY: rho_ocean, sa 889 878 890 USE averaging, &879 USE averaging, & 891 880 ONLY: rho_ocean_av, sa_av 892 881 893 USE control_parameters, &882 USE control_parameters, & 894 883 ONLY: average_count_3d 895 884 896 USE indices, &885 USE indices, & 897 886 ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt 898 887 … … 905 894 INTEGER(iwp) :: j !< loop index 906 895 INTEGER(iwp) :: k !< loop index 896 907 897 908 898 IF ( mode == 'allocate' ) THEN … … 936 926 DO j = nysg, nyng 937 927 DO k = nzb, nzt+1 938 rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) + & 939 rho_ocean(k,j,i) 928 rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) + rho_ocean(k,j,i) 940 929 ENDDO 941 930 ENDDO … … 964 953 965 954 CASE ( 'rho_sea_water' ) 966 IF ( ALLOCATED( rho_ocean_av ) ) THEN955 IF ( ALLOCATED( rho_ocean_av ) ) THEN 967 956 DO i = nxlg, nxrg 968 957 DO j = nysg, nyng 969 958 DO k = nzb, nzt+1 970 rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) / &959 rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) / & 971 960 REAL( average_count_3d, KIND=wp ) 972 961 ENDDO … … 976 965 977 966 CASE ( 'sa' ) 978 IF ( ALLOCATED( sa_av ) ) THEN967 IF ( ALLOCATED( sa_av ) ) THEN 979 968 DO i = nxlg, nxrg 980 969 DO j = nysg, nyng 981 970 DO k = nzb, nzt+1 982 sa_av(k,j,i) = sa_av(k,j,i) / & 983 REAL( average_count_3d, KIND=wp ) 971 sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d, KIND=wp ) 984 972 ENDDO 985 973 ENDDO … … 994 982 995 983 996 !------------------------------------------------------------------------------ !984 !--------------------------------------------------------------------------------------------------! 997 985 ! Description: 998 986 ! ------------ 999 987 !> Define 2D output variables. 1000 !------------------------------------------------------------------------------! 1001 SUBROUTINE ocean_data_output_2d( av, variable, found, grid, mode, local_pf, & 1002 nzb_do, nzt_do ) 1003 1004 USE arrays_3d, & 988 !--------------------------------------------------------------------------------------------------! 989 SUBROUTINE ocean_data_output_2d( av, variable, found, grid, mode, local_pf, nzb_do, nzt_do ) 990 991 USE arrays_3d, & 1005 992 ONLY: rho_ocean, sa 1006 993 1007 USE averaging, &994 USE averaging, & 1008 995 ONLY: rho_ocean_av, sa_av 1009 996 1010 USE indices, & 1011 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, & 1012 wall_flags_total_0 997 USE indices, & 998 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0 1013 999 1014 1000 IMPLICIT NONE … … 1031 1017 REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute 1032 1018 1033 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local 1034 !< array to which output data is resorted to 1019 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local array to which output data is resorted to 1035 1020 1036 1021 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable 1037 1022 1023 1038 1024 found = .TRUE. 1039 1025 resorted = .FALSE. … … 1048 1034 to_be_resorted => rho_ocean 1049 1035 ELSE 1050 IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN1036 IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN 1051 1037 ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1052 1038 rho_ocean_av = REAL( fill_value, KIND = wp ) … … 1059 1045 to_be_resorted => sa 1060 1046 ELSE 1061 IF ( .NOT. ALLOCATED( sa_av ) ) THEN1047 IF ( .NOT. ALLOCATED( sa_av ) ) THEN 1062 1048 ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1063 1049 sa_av = REAL( fill_value, KIND = wp ) … … 1077 1063 DO j = nys, nyn 1078 1064 DO k = nzb_do, nzt_do 1079 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 1080 REAL( fill_value, KIND = wp ), & 1081 BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 1065 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ), & 1066 BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 1082 1067 ENDDO 1083 1068 ENDDO … … 1085 1070 resorted = .TRUE. 1086 1071 ENDIF 1087 1072 1088 1073 END SUBROUTINE ocean_data_output_2d 1089 1074 1090 1091 !------------------------------------------------------------------------------ !1075 1076 !--------------------------------------------------------------------------------------------------! 1092 1077 ! Description: 1093 1078 ! ------------ 1094 1079 !> Define 3D output variables. 1095 !------------------------------------------------------------------------------ !1080 !--------------------------------------------------------------------------------------------------! 1096 1081 SUBROUTINE ocean_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 1097 1098 1099 USE arrays_3d, &1082 1083 1084 USE arrays_3d, & 1100 1085 ONLY: rho_ocean, sa 1101 1086 1102 USE averaging, &1087 USE averaging, & 1103 1088 ONLY: rho_ocean_av, sa_av 1104 1089 1105 USE indices, & 1106 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, & 1107 wall_flags_total_0 1090 USE indices, & 1091 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0 1108 1092 1109 1093 IMPLICIT NONE … … 1129 1113 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable 1130 1114 1115 1131 1116 found = .TRUE. 1132 1117 resorted = .FALSE. … … 1141 1126 to_be_resorted => rho_ocean 1142 1127 ELSE 1143 IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN1128 IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN 1144 1129 ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1145 1130 rho_ocean_av = REAL( fill_value, KIND = wp ) … … 1152 1137 to_be_resorted => sa 1153 1138 ELSE 1154 IF ( .NOT. ALLOCATED( sa_av ) ) THEN1139 IF ( .NOT. ALLOCATED( sa_av ) ) THEN 1155 1140 ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1156 1141 sa_av = REAL( fill_value, KIND = wp ) … … 1169 1154 DO j = nys, nyn 1170 1155 DO k = nzb_do, nzt_do 1171 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 1172 REAL( fill_value, KIND = wp ), & 1173 BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 1156 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ), & 1157 BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 1174 1158 ENDDO 1175 1159 ENDDO … … 1206 1190 1207 1191 1208 !------------------------------------------------------------------------------ !1192 !--------------------------------------------------------------------------------------------------! 1209 1193 ! Description: 1210 1194 ! ------------ 1211 1195 !> Header output for ocean parameters 1212 !------------------------------------------------------------------------------ !1196 !--------------------------------------------------------------------------------------------------! 1213 1197 SUBROUTINE ocean_header( io ) 1214 1198 1215 1216 1199 IMPLICIT NONE 1217 1200 1218 1201 INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file 1202 1219 1203 1220 1204 ! … … 1230 1214 ENDIF 1231 1215 1232 1 FORMAT (//' Ocean settings:'/ & 1233 ' ------------------------------------------'/) 1234 2 FORMAT (' --> Craik-Leibovich vortex force and Stokes drift switched', & 1235 ' on'/ & 1216 1 FORMAT (//' Ocean settings:'/' ------------------------------------------'/) 1217 2 FORMAT (' --> Craik-Leibovich vortex force and Stokes drift switched',' on'/ & 1236 1218 ' waveheight: ',F4.1,' m wavelength: ',F6.1,' m') 1237 3 FORMAT (' --> wave breaking generated turbulence switched on'/ & 1238 ' alpha: ',F4.1/ & 1239 ' timescale:',F5.1,' s') 1219 3 FORMAT (' --> wave breaking generated turbulence switched on'/ & 1220 ' alpha: ',F4.1/' timescale:',F5.1,' s') 1240 1221 4 FORMAT (' --> prognostic salinity equation is switched off' ) 1241 1222 5 FORMAT (' --> surface heat flux is switched off after ',F8.1,' s') … … 1244 1225 1245 1226 1246 !------------------------------------------------------------------------------ !1227 !--------------------------------------------------------------------------------------------------! 1247 1228 ! Description: 1248 1229 ! ------------ 1249 1230 !> Allocate arrays and assign pointers. 1250 !------------------------------------------------------------------------------ !1231 !--------------------------------------------------------------------------------------------------! 1251 1232 SUBROUTINE ocean_init_arrays 1252 1233 1253 USE indices, &1234 USE indices, & 1254 1235 ONLY: nxlg, nxrg, nyn, nyng, nys, nysg, nzb, nzt 1255 1236 1256 1237 IMPLICIT NONE 1257 1238 1258 ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &1259 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &1260 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &1239 ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 1240 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 1241 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 1261 1242 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1243 1262 1244 1263 1245 IF ( salinity ) THEN … … 1282 1264 1283 1265 prho => prho_1 1284 rho_ocean => rho_1 ! routines calc_mean_profile and diffusion_e require 1285 ! density to be a pointer 1266 rho_ocean => rho_1 ! routines calc_mean_profile and diffusion_e require density to be a pointer 1286 1267 1287 1268 ! … … 1296 1277 1297 1278 1298 !------------------------------------------------------------------------------ !1279 !--------------------------------------------------------------------------------------------------! 1299 1280 ! Description: 1300 1281 ! ------------ 1301 1282 !> Initialization of quantities needed for the ocean mode 1302 !------------------------------------------------------------------------------ !1283 !--------------------------------------------------------------------------------------------------! 1303 1284 SUBROUTINE ocean_init 1304 1285 1305 1286 1306 USE arrays_3d, &1307 ONLY: dzu, dzw, hyp, pt_init, ref_state, u_stokes_zu, u_stokes_zw, &1308 v_stokes_z u, v_stokes_zw, zu, zw1309 1310 USE basic_constants_and_equations_mod, &1287 USE arrays_3d, & 1288 ONLY: dzu, dzw, hyp, pt_init, ref_state, u_stokes_zu, u_stokes_zw, v_stokes_zu, & 1289 v_stokes_zw, zu, zw 1290 1291 USE basic_constants_and_equations_mod, & 1311 1292 ONLY: g 1312 1293 1313 USE basic_constants_and_equations_mod, &1294 USE basic_constants_and_equations_mod, & 1314 1295 ONLY: pi 1315 1296 1316 USE control_parameters, & 1317 ONLY: initializing_actions, molecular_viscosity, rho_surface, & 1318 rho_reference, surface_pressure, top_momentumflux_u, & 1319 top_momentumflux_v, use_single_reference_value 1320 1321 USE indices, & 1297 USE control_parameters, & 1298 ONLY: initializing_actions, molecular_viscosity, rho_surface, rho_reference, & 1299 surface_pressure, top_momentumflux_u, top_momentumflux_v, use_single_reference_value 1300 1301 USE indices, & 1322 1302 ONLY: nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt 1323 1303 1324 1304 USE kinds 1325 1305 1326 USE pegrid, &1306 USE pegrid, & 1327 1307 ONLY: myid 1328 1308 1329 USE statistics, &1309 USE statistics, & 1330 1310 ONLY: hom, statistic_regions 1331 1311 … … 1337 1317 INTEGER(iwp) :: n !< loop index 1338 1318 1339 REAL(wp) :: alpha !< angle of surface stress1340 REAL(wp) :: dum !< dummy argument1341 REAL(wp) :: pt_l !< local scalar for pt used in equation of state function1342 REAL(wp) :: sa_l !< local scalar for sa used in equation of state function1319 REAL(wp) :: alpha !< angle of surface stress 1320 REAL(wp) :: dum !< dummy argument 1321 REAL(wp) :: pt_l !< local scalar for pt used in equation of state function 1322 REAL(wp) :: sa_l !< local scalar for sa used in equation of state function 1343 1323 REAL(wp) :: velocity_amplitude !< local scalar for amplitude of Stokes drift velocity 1344 REAL(wp) :: x !< temporary variable to store surface stress along x1345 REAL(wp) :: y !< temporary variable to store surface stress along y1324 REAL(wp) :: x !< temporary variable to store surface stress along x 1325 REAL(wp) :: y !< temporary variable to store surface stress along y 1346 1326 1347 1327 REAL(wp), DIMENSION(nzb:nzt+1) :: rho_ocean_init !< local array for initial density … … 1351 1331 1352 1332 ! 1353 !-- In case of no restart run, calculate the inital salinity profile vcusing the1354 !-- g iven salinity gradients1333 !-- In case of no restart run, calculate the inital salinity profile using the given salinity 1334 !-- gradients. 1355 1335 IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN 1356 1336 1357 1337 sa_init = sa_surface 1358 1338 ! 1359 !-- Last arguments gives back the gradient at top level to be used as 1360 !-- possible Neumann boundary condition. This is not realized for the ocean 1361 !-- mode, therefore a dummy argument is used. 1339 !-- Last arguments gives back the gradient at top level to be used as possible Neumann boundary 1340 !-- condition. This is not realized for the ocean mode, therefore a dummy argument is used. 1362 1341 IF ( salinity ) THEN 1363 CALL init_vertical_profiles( sa_vertical_gradient_level_ind, &1364 sa_vertical_gradient_level, &1365 sa_vertical_gradient, sa_init, &1342 CALL init_vertical_profiles( sa_vertical_gradient_level_ind, & 1343 sa_vertical_gradient_level, & 1344 sa_vertical_gradient, sa_init, & 1366 1345 sa_surface, dum ) 1367 1346 ENDIF … … 1370 1349 ! 1371 1350 !-- Initialize required 3d-arrays 1372 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &1351 IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 1373 1352 TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN 1374 1353 ! … … 1395 1374 1396 1375 ! 1397 !-- Initialize new time levels (only done in order to set boundary values 1398 !-- including ghost points) 1376 !-- Initialize new time levels (only done in order to set boundary values including ghost points) 1399 1377 sa_p = sa 1400 1378 ! 1401 !-- Allthough tendency arrays are set in prognostic_equations, they have 1402 !-- have to be predefined here because they are used (but multiplied with 0) 1403 !-- there before they are set. 1379 !-- Allthough tendency arrays are set in prognostic_equations, they have to be predefined here 1380 !-- because they are used (but multiplied with 0) there before they are set. 1404 1381 tsa_m = 0.0_wp 1405 1382 … … 1416 1393 1417 1394 ! 1418 !-- Change sign of buoyancy/stability terms because density gradient is used 1419 !-- instead of the potential temperature gradient to calculate the buoyancy1395 !-- Change sign of buoyancy/stability terms because density gradient is used instead of the 1396 !-- potential temperature gradient to calculate the buoyancy. 1420 1397 atmos_ocean_sign = -1.0_wp 1421 1398 1422 1399 ! 1423 !-- Calculate initial vertical profile of hydrostatic pressure (in Pa) 1424 !-- and the reference density (used later in buoyancy term)1425 !-- First step: Calculate pressure using reference density 1400 !-- Calculate initial vertical profile of hydrostatic pressure (in Pa) and the reference density 1401 !-- (used later in buoyancy term). 1402 !-- First step: Calculate pressure using reference density. 1426 1403 hyp(nzt+1) = surface_pressure * 100.0_wp 1427 1404 hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1) … … 1435 1412 1436 1413 ! 1437 !-- Second step: Iteratively calculate in situ density (based on presssure) 1438 !-- and pressure(based on in situ density)1414 !-- Second step: Iteratively calculate in situ density (based on presssure) and pressure 1415 !-- (based on in situ density) 1439 1416 DO n = 1, 5 1440 1417 … … 1456 1433 hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1) 1457 1434 DO k = nzt-1, 0, -1 1458 hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k) & 1459 + rho_ocean_init(k+1) ) * dzu(k+1) 1435 hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k) + rho_ocean_init(k+1) ) * dzu(k+1) 1460 1436 ENDDO 1461 1437 … … 1470 1446 pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) ) 1471 1447 1472 prho_reference = prho_reference + dzu(k+1) * & 1473 eqn_state_seawater_func( 0.0_wp, pt_l, sa_l ) 1448 prho_reference = prho_reference + dzu(k+1) * eqn_state_seawater_func( 0.0_wp, pt_l, sa_l ) 1474 1449 1475 1450 ENDDO … … 1478 1453 1479 1454 ! 1480 !-- Calculate the 3d array of initial in situ and potential density, 1481 !-- based on the initial temperature and salinity profile1455 !-- Calculate the 3d array of initial in situ and potential density, based on the initial 1456 !-- temperature and salinity profile. 1482 1457 CALL eqn_state_seawater 1483 1458 … … 1520 1495 ENDIF 1521 1496 1522 velocity_amplitude = ( pi * stokes_waveheight / stokes_wavelength )**2 * &1497 velocity_amplitude = ( pi * stokes_waveheight / stokes_wavelength )**2 * & 1523 1498 SQRT( g * stokes_wavelength / ( 2.0_wp * pi ) ) 1524 1499 1525 1500 DO k = nzb, nzt 1526 u_stokes_zu(k) = velocity_amplitude * COS( alpha ) * &1501 u_stokes_zu(k) = velocity_amplitude * COS( alpha ) * & 1527 1502 EXP( 4.0_wp * pi * zu(k) / stokes_wavelength ) 1528 u_stokes_zw(k) = velocity_amplitude * COS( alpha ) * &1503 u_stokes_zw(k) = velocity_amplitude * COS( alpha ) * & 1529 1504 EXP( 4.0_wp * pi * zw(k) / stokes_wavelength ) 1530 v_stokes_zu(k) = velocity_amplitude * SIN( alpha ) * &1505 v_stokes_zu(k) = velocity_amplitude * SIN( alpha ) * & 1531 1506 EXP( 4.0_wp * pi * zu(k) / stokes_wavelength ) 1532 1507 v_stokes_zw(k) = velocity_amplitude * SIN( alpha ) * & … … 1545 1520 ! 1546 1521 !-- Calculate friction velocity at ocean surface 1547 u_star_wave_breaking = SQRT( SQRT( top_momentumflux_u**2 + & 1548 top_momentumflux_v**2 ) ) 1549 ! 1550 !-- Set the time scale of random forcing. The vertical grid spacing at the 1551 !-- ocean surface is assumed as the length scale of turbulence. 1522 u_star_wave_breaking = SQRT( SQRT( top_momentumflux_u**2 + top_momentumflux_v**2 ) ) 1523 ! 1524 !-- Set the time scale of random forcing. The vertical grid spacing at the ocean surface is 1525 !-- assumed as the length scale of turbulence. 1552 1526 !-- Formula follows Noh et al. (2004), JPO 1553 timescale_wave_breaking = 0.1_wp * dzw(nzt) / alpha_wave_breaking / & 1554 u_star_wave_breaking 1555 ! 1556 !-- Set random number seeds differently on the processor cores in order to 1557 !-- create different random number sequences 1527 timescale_wave_breaking = 0.1_wp * dzw(nzt) / alpha_wave_breaking / u_star_wave_breaking 1528 ! 1529 !-- Set random number seeds differently on the processor cores in order to create different 1530 !-- random number sequences 1558 1531 iran_ocean = iran_ocean + myid 1559 1532 ENDIF … … 1562 1535 1563 1536 1564 !------------------------------------------------------------------------------ !1537 !--------------------------------------------------------------------------------------------------! 1565 1538 ! Description: 1566 1539 ! ------------ 1567 1540 !> Call for all grid points 1568 !------------------------------------------------------------------------------ !1541 !--------------------------------------------------------------------------------------------------! 1569 1542 SUBROUTINE ocean_actions( location ) 1570 1543 … … 1588 1561 1589 1562 1590 !------------------------------------------------------------------------------ !1563 !--------------------------------------------------------------------------------------------------! 1591 1564 ! Description: 1592 1565 ! ------------ 1593 1566 !> Call for grid points i,j 1594 !------------------------------------------------------------------------------ !1567 !--------------------------------------------------------------------------------------------------! 1595 1568 SUBROUTINE ocean_actions_ij( i, j, location ) 1596 1569 1597 1598 INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction1599 INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction1600 1570 CHARACTER (LEN=*), INTENT(IN) :: location !< call location string 1601 INTEGER(iwp) :: dummy !< call location string 1571 1572 INTEGER(iwp) :: dummy !< call location string 1573 INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction 1574 INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction 1575 1602 1576 1603 1577 IF ( ocean_mode ) dummy = i + j … … 1619 1593 1620 1594 1621 !------------------------------------------------------------------------------ !1595 !--------------------------------------------------------------------------------------------------! 1622 1596 ! Description: 1623 1597 ! ------------ 1624 1598 !> Prognostic equation for salinity. 1625 1599 !> Vector-optimized version 1626 !------------------------------------------------------------------------------ !1600 !--------------------------------------------------------------------------------------------------! 1627 1601 SUBROUTINE ocean_prognostic_equations 1628 1602 1629 USE advec_s_bc_mod, &1603 USE advec_s_bc_mod, & 1630 1604 ONLY: advec_s_bc 1631 1605 1632 USE advec_s_pw_mod, &1606 USE advec_s_pw_mod, & 1633 1607 ONLY: advec_s_pw 1634 1608 1635 USE advec_s_up_mod, &1609 USE advec_s_up_mod, & 1636 1610 ONLY: advec_s_up 1637 1611 1638 USE advec_ws, &1612 USE advec_ws, & 1639 1613 ONLY: advec_s_ws 1640 1614 1641 USE arrays_3d, &1615 USE arrays_3d, & 1642 1616 ONLY: rdf_sc, tend, tsa_m 1643 1617 1644 USE control_parameters, & 1645 ONLY: bc_dirichlet_l, & 1646 bc_dirichlet_n, & 1647 bc_dirichlet_r, & 1648 bc_dirichlet_s, & 1649 bc_radiation_l, & 1650 bc_radiation_n, & 1651 bc_radiation_r, & 1652 bc_radiation_s, & 1653 dt_3d, intermediate_timestep_count, & 1654 intermediate_timestep_count_max, scalar_advec, simulated_time, & 1655 timestep_scheme, tsc, ws_scheme_sca 1656 1657 USE cpulog, & 1618 USE control_parameters, & 1619 ONLY: bc_dirichlet_l, & 1620 bc_dirichlet_n, & 1621 bc_dirichlet_r, & 1622 bc_dirichlet_s, & 1623 bc_radiation_l, & 1624 bc_radiation_n, & 1625 bc_radiation_r, & 1626 bc_radiation_s, & 1627 dt_3d, & 1628 intermediate_timestep_count, & 1629 intermediate_timestep_count_max, & 1630 scalar_advec, & 1631 simulated_time, & 1632 timestep_scheme, & 1633 tsc, & 1634 ws_scheme_sca 1635 1636 USE cpulog, & 1658 1637 ONLY: cpu_log, log_point_s 1659 1638 1660 USE diffusion_s_mod, &1639 USE diffusion_s_mod, & 1661 1640 ONLY: diffusion_s 1662 1641 … … 1669 1648 REAL(wp) :: sbt !< weighting factor for sub-time step 1670 1649 1650 1671 1651 ! 1672 1652 !-- Switch of the surface heat flux, if requested 1673 1653 IF ( surface_cooling_spinup_time /= 999999.9_wp ) THEN 1674 IF ( .NOT. surface_cooling_switched_off .AND. &1654 IF ( .NOT. surface_cooling_switched_off .AND. & 1675 1655 simulated_time >= surface_cooling_spinup_time ) THEN 1676 1656 … … 1682 1662 1683 1663 ! 1684 !-- Compute prognostic equations for the ocean mode 1685 !-- First, start with salinity 1664 !-- Compute prognostic equations for the ocean mode. 1665 !-- First, start with salinity. 1686 1666 IF ( salinity ) THEN 1687 1667 … … 1709 1689 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1710 1690 IF ( ws_scheme_sca ) THEN 1711 CALL advec_s_ws( advc_flags_s, sa, 'sa', &1712 bc_dirichlet_l .OR. bc_radiation_l, &1713 bc_dirichlet_n .OR. bc_radiation_n, &1714 bc_dirichlet_r .OR. bc_radiation_r, &1691 CALL advec_s_ws( advc_flags_s, sa, 'sa', & 1692 bc_dirichlet_l .OR. bc_radiation_l, & 1693 bc_dirichlet_n .OR. bc_radiation_n, & 1694 bc_dirichlet_r .OR. bc_radiation_r, & 1715 1695 bc_dirichlet_s .OR. bc_radiation_s ) 1716 1696 ELSE … … 1722 1702 ENDIF 1723 1703 1724 CALL diffusion_s( sa, &1725 surf_def_h(0)%sasws, surf_def_h(1)%sasws, &1726 surf_def_h(2)%sasws, &1727 surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws, &1728 surf_usm_h(0)%sasws, surf_usm_h(1)%sasws, &1729 surf_def_v(0)%sasws, surf_def_v(1)%sasws, &1730 surf_def_v(2)%sasws, surf_def_v(3)%sasws, &1731 surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, &1732 surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws, &1733 surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, &1704 CALL diffusion_s( sa, & 1705 surf_def_h(0)%sasws, surf_def_h(1)%sasws, & 1706 surf_def_h(2)%sasws, & 1707 surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws, & 1708 surf_usm_h(0)%sasws, surf_usm_h(1)%sasws, & 1709 surf_def_v(0)%sasws, surf_def_v(1)%sasws, & 1710 surf_def_v(2)%sasws, surf_def_v(3)%sasws, & 1711 surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, & 1712 surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws, & 1713 surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, & 1734 1714 surf_usm_v(2)%sasws, surf_usm_v(3)%sasws ) 1735 1715 … … 1743 1723 !DIR$ IVDEP 1744 1724 DO k = nzb+1, nzt 1745 sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & 1746 tsc(3) * tsa_m(k,j,i) ) & 1747 - tsc(5) * rdf_sc(k) * & 1748 ( sa(k,j,i) - sa_init(k) ) & 1749 ) & 1750 * MERGE( 1.0_wp, 0.0_wp, & 1751 BTEST( wall_flags_total_0(k,j,i), 0 ) & 1725 sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) ) & 1726 - tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) ) & 1727 ) & 1728 * MERGE( 1.0_wp, 0.0_wp, & 1729 BTEST( wall_flags_total_0(k,j,i), 0 ) & 1752 1730 ) 1753 1731 IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i) … … 1767 1745 ENDDO 1768 1746 ENDDO 1769 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) & 1770 THEN 1747 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 1771 1748 DO i = nxl, nxr 1772 1749 DO j = nys, nyn 1773 1750 DO k = nzb+1, nzt 1774 tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & 1775 5.3125_wp * tsa_m(k,j,i) 1751 tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tsa_m(k,j,i) 1776 1752 ENDDO 1777 1753 ENDDO … … 1793 1769 1794 1770 1795 !------------------------------------------------------------------------------ !1771 !--------------------------------------------------------------------------------------------------! 1796 1772 ! Description: 1797 1773 ! ------------ 1798 1774 !> Prognostic equations for ocean mode (so far, salinity only) 1799 1775 !> Cache-optimized version 1800 !------------------------------------------------------------------------------ !1776 !--------------------------------------------------------------------------------------------------! 1801 1777 SUBROUTINE ocean_prognostic_equations_ij( i, j, i_omp_start, tn ) 1802 1778 1803 USE advec_s_pw_mod, &1779 USE advec_s_pw_mod, & 1804 1780 ONLY: advec_s_pw 1805 1781 1806 USE advec_s_up_mod, &1782 USE advec_s_up_mod, & 1807 1783 ONLY: advec_s_up 1808 1784 1809 USE advec_ws, &1785 USE advec_ws, & 1810 1786 ONLY: advec_s_ws 1811 1787 1812 USE arrays_3d, &1788 USE arrays_3d, & 1813 1789 ONLY: diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, rdf_sc, tend, tsa_m 1814 1790 1815 USE control_parameters, & 1816 ONLY: bc_dirichlet_l, & 1817 bc_dirichlet_n, & 1818 bc_dirichlet_r, & 1819 bc_dirichlet_s, & 1820 bc_radiation_l, & 1821 bc_radiation_n, & 1822 bc_radiation_r, & 1823 bc_radiation_s, & 1824 dt_3d, intermediate_timestep_count, & 1825 intermediate_timestep_count_max, simulated_time, & 1826 timestep_scheme, tsc, ws_scheme_sca 1827 1828 USE diffusion_s_mod, & 1791 USE control_parameters, & 1792 ONLY: bc_dirichlet_l, & 1793 bc_dirichlet_n, & 1794 bc_dirichlet_r, & 1795 bc_dirichlet_s, & 1796 bc_radiation_l, & 1797 bc_radiation_n, & 1798 bc_radiation_r, & 1799 bc_radiation_s, & 1800 dt_3d, & 1801 intermediate_timestep_count, & 1802 intermediate_timestep_count_max, & 1803 simulated_time, & 1804 timestep_scheme, & 1805 tsc, & 1806 ws_scheme_sca 1807 1808 USE diffusion_s_mod, & 1829 1809 ONLY: diffusion_s 1830 1810 … … 1832 1812 1833 1813 INTEGER(iwp) :: i !< loop index x direction 1834 INTEGER(iwp) :: i_omp_start !< first loop index of i-loop in calling & 1835 !< routine prognostic_equations 1814 INTEGER(iwp) :: i_omp_start !< first loop index of i-loop in calling routine prognostic_equations 1836 1815 INTEGER(iwp) :: j !< loop index y direction 1837 1816 INTEGER(iwp) :: k !< loop index z direction … … 1842 1821 !-- Switch of the surface heat flux, if requested 1843 1822 IF ( surface_cooling_spinup_time /= 999999.9_wp ) THEN 1844 IF ( .NOT. surface_cooling_switched_off .AND. &1823 IF ( .NOT. surface_cooling_switched_off .AND. & 1845 1824 simulated_time >= surface_cooling_spinup_time ) THEN 1846 1825 … … 1852 1831 1853 1832 ! 1854 !-- Compute prognostic equations for the ocean mode 1855 !-- First, start with tendency-terms for salinity 1833 !-- Compute prognostic equations for the ocean mode. 1834 !-- First, start with tendency-terms for salinity. 1856 1835 IF ( salinity ) THEN 1857 1836 1858 1837 tend(:,j,i) = 0.0_wp 1859 IF ( timestep_scheme(1:5) == 'runge' ) & 1860 THEN 1838 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1861 1839 IF ( ws_scheme_sca ) THEN 1862 CALL advec_s_ws( advc_flags_s, &1863 i, j, sa, 'sa', flux_s_sa, diss_s_sa, flux_l_sa, &1864 diss_l_sa, i_omp_start, tn,&1865 bc_dirichlet_l .OR. bc_radiation_l, &1866 bc_dirichlet_n .OR. bc_radiation_n, &1867 bc_dirichlet_r .OR. bc_radiation_r, &1840 CALL advec_s_ws( advc_flags_s, & 1841 i, j, sa, 'sa', flux_s_sa, diss_s_sa, flux_l_sa, diss_l_sa, & 1842 i_omp_start, tn, & 1843 bc_dirichlet_l .OR. bc_radiation_l, & 1844 bc_dirichlet_n .OR. bc_radiation_n, & 1845 bc_dirichlet_r .OR. bc_radiation_r, & 1868 1846 bc_dirichlet_s .OR. bc_radiation_s ) 1869 1847 ELSE … … 1873 1851 CALL advec_s_up( i, j, sa ) 1874 1852 ENDIF 1875 CALL diffusion_s( i, j, sa, & 1876 surf_def_h(0)%sasws, surf_def_h(1)%sasws, & 1877 surf_def_h(2)%sasws, & 1878 surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws, & 1879 surf_usm_h(0)%sasws, surf_usm_h(1)%sasws, & 1880 surf_def_v(0)%sasws, surf_def_v(1)%sasws, & 1881 surf_def_v(2)%sasws, surf_def_v(3)%sasws, & 1882 surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, & 1883 surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws, & 1884 surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, & 1885 surf_usm_v(2)%sasws, surf_usm_v(3)%sasws ) 1853 CALL diffusion_s( i, j, sa, & 1854 surf_def_h(0)%sasws, surf_def_h(1)%sasws, surf_def_h(2)%sasws, & 1855 surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws, & 1856 surf_usm_h(0)%sasws, surf_usm_h(1)%sasws, & 1857 surf_def_v(0)%sasws, surf_def_v(1)%sasws, surf_def_v(2)%sasws, & 1858 surf_def_v(3)%sasws, & 1859 surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, surf_lsm_v(2)%sasws, & 1860 surf_lsm_v(3)%sasws, & 1861 surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, surf_usm_v(2)%sasws, & 1862 surf_usm_v(3)%sasws ) 1886 1863 1887 1864 ! CALL user_actions( i, j, 'sa-tendency' ) ToDo: find general solution for dependency between modules … … 1891 1868 DO k = nzb+1, nzt 1892 1869 1893 sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * & 1894 ( tsc(2) * tend(k,j,i) + & 1895 tsc(3) * tsa_m(k,j,i) ) & 1896 - tsc(5) * rdf_sc(k) & 1897 * ( sa(k,j,i) - sa_init(k) ) & 1898 ) * MERGE( 1.0_wp, 0.0_wp, & 1899 BTEST( wall_flags_total_0(k,j,i), 0 ) ) 1870 sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) ) & 1871 - tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) ) & 1872 ) * MERGE( 1.0_wp, 0.0_wp, & 1873 BTEST( wall_flags_total_0(k,j,i), 0 ) & 1874 ) 1900 1875 1901 1876 IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i) … … 1910 1885 tsa_m(k,j,i) = tend(k,j,i) 1911 1886 ENDDO 1912 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) & 1913 THEN 1887 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 1914 1888 DO k = nzb+1, nzt 1915 tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & 1916 5.3125_wp * tsa_m(k,j,i) 1889 tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tsa_m(k,j,i) 1917 1890 ENDDO 1918 1891 ENDIF … … 1927 1900 END SUBROUTINE ocean_prognostic_equations_ij 1928 1901 1929 !------------------------------------------------------------------------------ !1902 !--------------------------------------------------------------------------------------------------! 1930 1903 ! Description: 1931 1904 ! ------------ 1932 1905 !> Boundary conditions for ocean model 1933 !------------------------------------------------------------------------------ !1906 !--------------------------------------------------------------------------------------------------! 1934 1907 SUBROUTINE ocean_boundary_conditions 1935 1908 … … 1943 1916 1944 1917 ! 1945 !-- Boundary conditions for salinity 1946 !-- Bottom boundary: Neumann condition because salinity flux is always 1947 !-- given. 1918 !-- Boundary conditions for salinity. 1919 !-- Bottom boundary: Neumann condition because salinity flux is always given. 1948 1920 DO l = 0, 1 1949 1921 !$OMP PARALLEL DO PRIVATE( i, j, k ) … … 1965 1937 END SUBROUTINE ocean_boundary_conditions 1966 1938 1967 !------------------------------------------------------------------------------ !1939 !--------------------------------------------------------------------------------------------------! 1968 1940 ! Description: 1969 1941 ! ------------ 1970 1942 !> Swapping of timelevels. 1971 !------------------------------------------------------------------------------ !1943 !--------------------------------------------------------------------------------------------------! 1972 1944 SUBROUTINE ocean_swap_timelevel( mod_count ) 1973 1945 … … 1994 1966 1995 1967 1996 !------------------------------------------------------------------------------ !1968 !--------------------------------------------------------------------------------------------------! 1997 1969 ! Description: 1998 1970 ! ------------ 1999 1971 !> Read module-specific global restart data (Fortran binary format). 2000 !------------------------------------------------------------------------------ !1972 !--------------------------------------------------------------------------------------------------! 2001 1973 SUBROUTINE ocean_rrd_global_ftn( found ) 2002 1974 2003 1975 2004 USE control_parameters, &1976 USE control_parameters, & 2005 1977 ONLY: length, restart_string 2006 1978 … … 2063 2035 2064 2036 2065 !------------------------------------------------------------------------------ !2037 !--------------------------------------------------------------------------------------------------! 2066 2038 ! Description: 2067 2039 ! ------------ 2068 2040 !> Read module-specific global restart data (MPI-IO). 2069 !------------------------------------------------------------------------------ !2041 !--------------------------------------------------------------------------------------------------! 2070 2042 SUBROUTINE ocean_rrd_global_mpi 2071 2043 … … 2088 2060 2089 2061 2090 !------------------------------------------------------------------------------ !2062 !--------------------------------------------------------------------------------------------------! 2091 2063 ! Description: 2092 2064 ! ------------ 2093 2065 !> Read module-specific local restart data arrays (Fortran binary format). 2094 !------------------------------------------------------------------------------! 2095 SUBROUTINE ocean_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, & 2096 nxr_on_file, nynf, nync, nyn_on_file, nysf, & 2097 nysc, nys_on_file, tmp_3d, found ) 2098 2099 USE averaging, & 2066 !--------------------------------------------------------------------------------------------------! 2067 SUBROUTINE ocean_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync, & 2068 nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found ) 2069 2070 USE averaging, & 2100 2071 ONLY: rho_ocean_av, sa_av 2101 2072 2102 USE control_parameters, &2073 USE control_parameters, & 2103 2074 ONLY: length, restart_string 2104 2075 2105 USE indices, &2076 USE indices, & 2106 2077 ONLY: nbgp, nxlg, nxrg, nyng, nysg, nzb, nzt 2107 2078 … … 2127 2098 LOGICAL, INTENT(OUT) :: found 2128 2099 2129 REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !< 2100 REAL(wp), & 2101 DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) & 2102 :: tmp_3d !< 2130 2103 2131 2104 … … 2139 2112 ENDIF 2140 2113 IF ( k == 1 ) READ ( 13 ) tmp_3d 2141 rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &2142 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)2114 rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 2115 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 2143 2116 2144 2117 CASE ( 'sa' ) 2145 2118 IF ( k == 1 ) READ ( 13 ) tmp_3d 2146 sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &2147 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)2119 sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 2120 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 2148 2121 2149 2122 CASE ( 'sa_av' ) … … 2152 2125 ENDIF 2153 2126 IF ( k == 1 ) READ ( 13 ) tmp_3d 2154 sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &2155 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)2127 sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 2128 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 2156 2129 2157 2130 CASE DEFAULT … … 2163 2136 2164 2137 2165 !------------------------------------------------------------------------------ !2138 !--------------------------------------------------------------------------------------------------! 2166 2139 ! Description: 2167 2140 ! ------------ 2168 2141 !> Read module-specific local restart data arrays (MPI-IO). 2169 !------------------------------------------------------------------------------ !2142 !--------------------------------------------------------------------------------------------------! 2170 2143 SUBROUTINE ocean_rrd_local_mpi 2171 2144 2172 USE averaging, &2145 USE averaging, & 2173 2146 ONLY: rho_ocean_av, sa_av 2174 2147 2175 USE indices, &2148 USE indices, & 2176 2149 ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt 2177 2150 … … 2183 2156 CALL rd_mpi_io_check_array( 'rho_ocean_av' , found = array_found ) 2184 2157 IF ( array_found ) THEN 2185 IF ( .NOT. ALLOCATED( rho_ocean_av ) ) ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 2158 IF ( .NOT. ALLOCATED( rho_ocean_av ) ) & 2159 ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 2186 2160 CALL rrd_mpi_io( 'rho_ocean_av', rho_ocean_av ) 2187 2161 ENDIF … … 2198 2172 2199 2173 2200 !------------------------------------------------------------------------------ !2174 !--------------------------------------------------------------------------------------------------! 2201 2175 ! Description: 2202 2176 ! ------------ 2203 2177 !> This routine writes the respective restart data for the ocean module. 2204 !------------------------------------------------------------------------------ !2178 !--------------------------------------------------------------------------------------------------! 2205 2179 SUBROUTINE ocean_wrd_global 2206 2180 … … 2258 2232 CALL wrd_mpi_io_global_array( 'sa_vertical_gradient', sa_vertical_gradient ) 2259 2233 CALL wrd_mpi_io_global_array( 'sa_vertical_gradient_level', sa_vertical_gradient_level ) 2260 CALL wrd_mpi_io_global_array( 'sa_vertical_gradient_level_ind', sa_vertical_gradient_level_ind ) 2234 CALL wrd_mpi_io_global_array( 'sa_vertical_gradient_level_ind', & 2235 sa_vertical_gradient_level_ind ) 2261 2236 CALL wrd_mpi_io( 'stokes_waveheight', stokes_waveheight ) 2262 2237 CALL wrd_mpi_io( 'stokes_wavelength', stokes_wavelength ) … … 2271 2246 2272 2247 2273 !------------------------------------------------------------------------------ !2248 !--------------------------------------------------------------------------------------------------! 2274 2249 ! Description: 2275 2250 ! ------------ 2276 2251 !> This routine writes the respective restart data for the ocean module. 2277 !------------------------------------------------------------------------------ !2252 !--------------------------------------------------------------------------------------------------! 2278 2253 SUBROUTINE ocean_wrd_local 2279 2254 2280 USE averaging, &2255 USE averaging, & 2281 2256 ONLY: rho_ocean_av, sa_av 2282 2257 … … 2310 2285 2311 2286 2312 !------------------------------------------------------------------------------ !2313 ! Description: 2314 ! ------------ 2315 !> This routine calculates the Craik Leibovich vortex force and the additional 2316 !> effect of the Stokes drift on the Coriolis force2287 !--------------------------------------------------------------------------------------------------! 2288 ! Description: 2289 ! ------------ 2290 !> This routine calculates the Craik Leibovich vortex force and the additional effect of the Stokes 2291 !> drift on the Coriolis force. 2317 2292 !> Call for all gridpoints. 2318 !------------------------------------------------------------------------------ !2293 !--------------------------------------------------------------------------------------------------! 2319 2294 SUBROUTINE stokes_drift_terms( component ) 2320 2295 2321 USE arrays_3d, & 2322 ONLY: ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, & 2323 v_stokes_zw, w, tend 2324 2325 USE basic_constants_and_equations_mod, & 2296 USE arrays_3d, & 2297 ONLY: ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, v_stokes_zw, w, tend 2298 2299 USE basic_constants_and_equations_mod, & 2326 2300 ONLY: pi 2327 2301 2328 USE control_parameters, &2302 USE control_parameters, & 2329 2303 ONLY: f, fs, message_string, rotation_angle 2330 2304 2331 USE grid_variables, &2305 USE grid_variables, & 2332 2306 ONLY: ddx, ddy 2333 2307 2334 USE indices, &2308 USE indices, & 2335 2309 ONLY: nxl, nxr, nys, nysv, nyn, nzb, nzt 2336 2310 … … 2344 2318 REAL(wp) :: cos_rot_angle !< cosine of model rotation angle 2345 2319 REAL(wp) :: sin_rot_angle !< sine of model rotation angle 2320 2346 2321 2347 2322 ! … … 2355 2330 DO j = nysv, nyn 2356 2331 DO k = nzb+1, nzt 2357 tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * ( &2358 0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)&2359 + v(k,j,i) - v(k,j,i-1) ) * ddx&2360 - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) ) * ddy&2361 ) &2332 tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * ( & 2333 0.5 * ( v(k,j+1,i) - v(k,j+1,i-1) & 2334 + v(k,j,i) - v(k,j,i-1) ) * ddx & 2335 - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) ) * ddy & 2336 ) & 2362 2337 + f * v_stokes_zu(k) 2363 2338 ENDDO … … 2371 2346 DO j = nysv, nyn 2372 2347 DO k = nzb+1, nzt 2373 tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * ( &2374 0.5 * ( v(k,j,i+1) - v(k,j,i-1) ) * ddx&2375 - 0.5 * ( u(k,j,i) - u(k,j-1,i)&2376 + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy&2348 tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * ( & 2349 0.5 * ( v(k,j,i+1) - v(k,j,i-1) ) * ddx & 2350 - 0.5 * ( u(k,j,i) - u(k,j-1,i) & 2351 + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 2377 2352 ) & 2378 2353 - f * u_stokes_zu(k) … … 2393 2368 DO j = nys, nyn 2394 2369 DO k = nzb+1, nzt 2395 tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * ( & 2396 0.5 * ( u(k+1,j,i) - u(k,j,i) & 2397 + u(k+1,j,i+1) - u(k,j,i+1) & 2398 ) * ddzu(k+1) & 2399 - 0.5 * ( w(k,j,i+1) - w(k,j,i-1) & 2400 ) * ddx ) & 2401 - v_stokes_zw(k) * ( & 2402 0.5 * ( w(k,j+1,i) - w(k,j-1,i) & 2403 ) * ddy & 2404 - 0.5 * ( v(k+1,j,i) - v(k,j,i) & 2405 + v(k+1,j+1,i) - v(k,j+1,i) & 2406 ) * ddzu(k) ) & 2407 + fs * ( & 2408 sin_rot_angle * v_stokes_zw(k) & 2409 + cos_rot_angle * u_stokes_zw(k) & 2410 ) 2370 tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * ( & 2371 0.5 * ( u(k+1,j,i) - u(k,j,i) & 2372 + u(k+1,j,i+1) - u(k,j,i+1) & 2373 ) * ddzu(k+1) & 2374 - 0.5 * ( w(k,j,i+1) - w(k,j,i-1) & 2375 ) * ddx ) & 2376 - v_stokes_zw(k) * ( & 2377 0.5 * ( w(k,j+1,i) - w(k,j-1,i) & 2378 ) * ddy & 2379 - 0.5 * ( v(k+1,j,i) - v(k,j,i) & 2380 + v(k+1,j+1,i) - v(k,j+1,i) & 2381 ) * ddzu(k) ) & 2382 + fs * ( sin_rot_angle * v_stokes_zw(k) & 2383 + cos_rot_angle * u_stokes_zw(k) & 2384 ) 2411 2385 ENDDO 2412 2386 ENDDO … … 2414 2388 2415 2389 CASE DEFAULT 2416 WRITE( message_string, * ) 'wrong component of Stokes force: ', & 2417 component 2390 WRITE( message_string, * ) 'wrong component of Stokes force: ', component 2418 2391 CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 ) 2419 2392 … … 2423 2396 2424 2397 2425 !------------------------------------------------------------------------------ !2426 ! Description: 2427 ! ------------ 2428 !> This routine calculates the Craik Leibovich vortex force and the additional 2429 !> effect of the Stokes drift on the Coriolis force2398 !--------------------------------------------------------------------------------------------------! 2399 ! Description: 2400 ! ------------ 2401 !> This routine calculates the Craik Leibovich vortex force and the additional effect of the Stokes 2402 !> drift on the Coriolis force. 2430 2403 !> Call for gridpoints i,j. 2431 !------------------------------------------------------------------------------ !2404 !--------------------------------------------------------------------------------------------------! 2432 2405 2433 2406 SUBROUTINE stokes_drift_terms_ij( i, j, component ) 2434 2407 2435 USE arrays_3d, & 2436 ONLY: ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, & 2437 v_stokes_zw, w, tend 2438 2439 USE basic_constants_and_equations_mod, & 2408 USE arrays_3d, & 2409 ONLY: ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, v_stokes_zw, w, tend 2410 2411 USE basic_constants_and_equations_mod, & 2440 2412 ONLY: pi 2441 2413 2442 USE control_parameters, &2414 USE control_parameters, & 2443 2415 ONLY: f, fs, message_string, rotation_angle 2444 2416 2445 USE grid_variables, &2417 USE grid_variables, & 2446 2418 ONLY: ddx, ddy 2447 2419 2448 USE indices, &2420 USE indices, & 2449 2421 ONLY: nzb, nzt 2450 2422 … … 2459 2431 REAL(wp) :: sin_rot_angle !< sine of model rotation angle 2460 2432 2433 2461 2434 ! 2462 2435 !-- Compute Stokes terms for the respective velocity components … … 2467 2440 CASE ( 1 ) 2468 2441 DO k = nzb+1, nzt 2469 tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * ( &2470 0.5 * ( v(k,j+1,i) - v(k,j+1,i-1)&2471 + v(k,j,i) - v(k,j,i-1) ) * ddx&2472 - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) ) * ddy&2473 ) &2442 tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * ( & 2443 0.5 * ( v(k,j+1,i) - v(k,j+1,i-1) & 2444 + v(k,j,i) - v(k,j,i-1) ) * ddx & 2445 - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) ) * ddy & 2446 ) & 2474 2447 + f * v_stokes_zu(k) 2475 2448 ENDDO … … 2478 2451 CASE ( 2 ) 2479 2452 DO k = nzb+1, nzt 2480 tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * ( &2481 0.5 * ( v(k,j,i+1) - v(k,j,i-1) ) * ddx&2482 - 0.5 * ( u(k,j,i) - u(k,j-1,i)&2483 + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy&2484 ) &2453 tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * ( & 2454 0.5 * ( v(k,j,i+1) - v(k,j,i-1) ) * ddx & 2455 - 0.5 * ( u(k,j,i) - u(k,j-1,i) & 2456 + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 2457 ) & 2485 2458 - f * u_stokes_zu(k) 2486 2459 ENDDO … … 2496 2469 2497 2470 DO k = nzb+1, nzt 2498 tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * ( & 2499 0.5 * ( u(k+1,j,i) - u(k,j,i) & 2500 + u(k+1,j,i+1) - u(k,j,i+1) & 2501 ) * ddzu(k+1) & 2502 - 0.5 * ( w(k,j,i+1) - w(k,j,i-1) & 2503 ) * ddx ) & 2504 - v_stokes_zw(k) * ( & 2505 0.5 * ( w(k,j+1,i) - w(k,j-1,i) & 2506 ) * ddy & 2507 - 0.5 * ( v(k+1,j,i) - v(k,j,i) & 2508 + v(k+1,j+1,i) - v(k,j+1,i) & 2509 ) * ddzu(k) ) & 2510 + fs * ( sin_rot_angle * v_stokes_zw(k) & 2511 + cos_rot_angle * u_stokes_zw(k) & 2471 tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * ( 0.5 * ( u(k+1,j,i) - u(k,j,i) & 2472 + u(k+1,j,i+1) - u(k,j,i+1) & 2473 ) * ddzu(k+1) & 2474 - 0.5 * ( w(k,j,i+1) - w(k,j,i-1) & 2475 ) * ddx ) & 2476 - v_stokes_zw(k) * ( 0.5 * ( w(k,j+1,i) - w(k,j-1,i) & 2477 ) * ddy & 2478 - 0.5 * ( v(k+1,j,i) - v(k,j,i) & 2479 + v(k+1,j+1,i) - v(k,j+1,i) & 2480 ) * ddzu(k) ) & 2481 + fs * ( sin_rot_angle * v_stokes_zw(k) & 2482 + cos_rot_angle * u_stokes_zw(k) & 2512 2483 ) 2513 2484 ENDDO … … 2522 2493 2523 2494 2524 !------------------------------------------------------------------------------ !2525 ! Description: 2526 ! ------------ 2527 !> This routine calculates turbulence generated by wave breaking near the ocean 2528 !> surface, followinga parameterization given in Noh et al. (2004), JPO2495 !--------------------------------------------------------------------------------------------------! 2496 ! Description: 2497 ! ------------ 2498 !> This routine calculates turbulence generated by wave breaking near the ocean surface, following 2499 !> a parameterization given in Noh et al. (2004), JPO 2529 2500 !> Call for all gridpoints. 2530 !> TODO: so far, this routine only works if the model time step has about the 2531 !> s ame value as the time scale of wave breaking!2532 !------------------------------------------------------------------------------ !2501 !> TODO: so far, this routine only works if the model time step has about the same value as the time 2502 !> scale of wave breaking! 2503 !--------------------------------------------------------------------------------------------------! 2533 2504 SUBROUTINE wave_breaking_term( component ) 2534 2505 2535 USE arrays_3d, &2506 USE arrays_3d, & 2536 2507 ONLY: u_p, v_p 2537 2508 2538 USE control_parameters, &2509 USE control_parameters, & 2539 2510 ONLY: dt_3d, message_string 2540 2511 2541 USE indices, &2512 USE indices, & 2542 2513 ONLY: nxl, nxlu, nxr, nys, nysv, nyn, nzt 2543 2514 … … 2548 2519 INTEGER(iwp) :: j !< loop index along y 2549 2520 2550 REAL(wp) :: random_gauss !< function that creates a random number with a 2551 !< Gaussian distribution 2521 REAL(wp) :: random_gauss !< function that creates a random number with a Gaussian distribution 2552 2522 2553 2523 … … 2562 2532 DO i = nxlu, nxr 2563 2533 DO j = nys, nyn 2564 u_p(nzt,j,i) = u_p(nzt,j,i) + &2565 ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &2566 * alpha_wave_breaking * u_star_wave_breaking &2534 u_p(nzt,j,i) = u_p(nzt,j,i) + & 2535 ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) & 2536 * alpha_wave_breaking * u_star_wave_breaking & 2567 2537 / timescale_wave_breaking * dt_3d 2568 2538 ENDDO … … 2573 2543 DO i = nxl, nxr 2574 2544 DO j = nysv, nyn 2575 v_p(nzt,j,i) = v_p(nzt,j,i) + &2576 ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &2577 * alpha_wave_breaking * u_star_wave_breaking &2545 v_p(nzt,j,i) = v_p(nzt,j,i) + & 2546 ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) & 2547 * alpha_wave_breaking * u_star_wave_breaking & 2578 2548 / timescale_wave_breaking * dt_3d 2579 2549 ENDDO … … 2581 2551 2582 2552 CASE DEFAULT 2583 WRITE( message_string, * ) 'wrong component of wave breaking: ', & 2584 component 2553 WRITE( message_string, * ) 'wrong component of wave breaking: ', component 2585 2554 CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 ) 2586 2555 … … 2590 2559 2591 2560 2592 !------------------------------------------------------------------------------ !2593 ! Description: 2594 ! ------------ 2595 !> This routine calculates turbulence generated by wave breaking near the ocean 2596 !> surface, following aparameterization given in Noh et al. (2004), JPO2561 !--------------------------------------------------------------------------------------------------! 2562 ! Description: 2563 ! ------------ 2564 !> This routine calculates turbulence generated by wave breaking near the ocean surface, following a 2565 !> parameterization given in Noh et al. (2004), JPO 2597 2566 !> Call for gridpoint i,j. 2598 !> TODO: so far, this routine only works if the model time step has about the 2599 !> same value as the timescale of wave breaking!2600 !------------------------------------------------------------------------------ !2567 !> TODO: so far, this routine only works if the model time step has about the same value as the time 2568 !> scale of wave breaking! 2569 !--------------------------------------------------------------------------------------------------! 2601 2570 SUBROUTINE wave_breaking_term_ij( i, j, component ) 2602 2571 2603 USE arrays_3d, &2572 USE arrays_3d, & 2604 2573 ONLY: u_p, v_p 2605 2574 2606 USE control_parameters, &2575 USE control_parameters, & 2607 2576 ONLY: dt_3d, message_string 2608 2577 2609 USE indices, &2578 USE indices, & 2610 2579 ONLY: nzt 2611 2580 … … 2616 2585 INTEGER(iwp) :: j !< loop index along y 2617 2586 2618 REAL(wp) :: random_gauss !< function that creates a random number with a 2619 !< Gaussian distribution 2587 REAL(wp) :: random_gauss !< function that creates a random number with a Gaussian distribution 2620 2588 2621 2589 ! … … 2626 2594 !-- u-/v-component 2627 2595 CASE ( 1 ) 2628 u_p(nzt,j,i) = u_p(nzt,j,i) + &2629 ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &2630 * alpha_wave_breaking * u_star_wave_breaking &2596 u_p(nzt,j,i) = u_p(nzt,j,i) + & 2597 ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) & 2598 * alpha_wave_breaking * u_star_wave_breaking & 2631 2599 / timescale_wave_breaking * dt_3d 2632 2600 2633 2601 CASE ( 2 ) 2634 v_p(nzt,j,i) = v_p(nzt,j,i) + &2635 ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) &2636 * alpha_wave_breaking * u_star_wave_breaking &2602 v_p(nzt,j,i) = v_p(nzt,j,i) + & 2603 ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) & 2604 * alpha_wave_breaking * u_star_wave_breaking & 2637 2605 / timescale_wave_breaking * dt_3d 2638 2606 2639 2607 CASE DEFAULT 2640 WRITE( message_string, * ) 'wrong component of wave breaking: ', & 2641 component 2608 WRITE( message_string, * ) 'wrong component of wave breaking: ', component 2642 2609 CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 ) 2643 2610
Note: See TracChangeset
for help on using the changeset viewer.
