!> @file ocean_mod.f90
!--------------------------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
! Public License as published by the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
! Public License for more details.
!
! You should have received a copy of the GNU General Public License along with PALM. If not, see
! <http://www.gnu.org/licenses/>.
!
! Copyright 2017-2021 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------------------------!
!
! Authors:
! --------
! @author Siegfried Raasch
!
! Description:
! ------------
!> Inititalizes the vertical profiles of scalar quantities.
!--------------------------------------------------------------------------------------------------!
 SUBROUTINE init_vertical_profiles( vertical_gradient_level_ind,                                   &
                                    vertical_gradient_level,                                       &
                                    vertical_gradient, initial_profile,                            &
                                    surface_value, bc_top_gradient )

    USE arrays_3d,                                                                                 &
        ONLY:  dzu, zu

    USE control_parameters,                                                                        &
        ONLY:  ocean_mode

    USE indices,                                                                                   &
        ONLY:  nz, nzt

    USE kinds

    IMPLICIT NONE

    INTEGER(iwp) ::  i  !< loop counter
    INTEGER(iwp) ::  k  !< loop counter

    INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind  !< vertical grid indices for gradient levels

    REAL(wp)     ::  bc_top_gradient  !< model top gradient
    REAL(wp)     ::  gradient         !< vertica gradient of the respective quantity
    REAL(wp)     ::  surface_value    !< surface value of the respecitve quantity

    REAL(wp), DIMENSION(0:nz+1) ::  initial_profile          !< initialisation profile
    REAL(wp), DIMENSION(1:10)   ::  vertical_gradient        !< given vertical gradient
    REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level  !< given vertical gradient level

    i = 1
    gradient = 0.0_wp

    IF ( .NOT. ocean_mode )  THEN

       vertical_gradient_level_ind(1) = 0
       DO  k = 1, nzt+1
          IF ( i < 11 )  THEN
             IF ( vertical_gradient_level(i) < zu(k)  .AND.                                        &
                  vertical_gradient_level(i) >= 0.0_wp )  THEN
                gradient = vertical_gradient(i) / 100.0_wp
                vertical_gradient_level_ind(i) = k - 1
                i = i + 1
             ENDIF
          ENDIF
          IF ( gradient /= 0.0_wp )  THEN
             IF ( k /= 1 )  THEN
                initial_profile(k) = initial_profile(k-1) + dzu(k) * gradient
             ELSE
                initial_profile(k) = initial_profile(k-1) + dzu(k) * gradient
             ENDIF
          ELSE
             initial_profile(k) = initial_profile(k-1)
          ENDIF
!
!--       Avoid negative values of scalars
          IF ( initial_profile(k) < 0.0_wp )  THEN
             initial_profile(k) = 0.0_wp
          ENDIF
       ENDDO

    ELSE

!
!--    In ocean mode, profiles are constructed starting from the ocean surface, which is at the top
!--    of the model domain
       vertical_gradient_level_ind(1) = nzt+1
       DO  k = nzt, 0, -1
          IF ( i < 11 )  THEN
             IF ( vertical_gradient_level(i) > zu(k)  .AND.                                        &
                  vertical_gradient_level(i) <= 0.0_wp )  THEN
                gradient = vertical_gradient(i) / 100.0_wp
                vertical_gradient_level_ind(i) = k + 1
                i = i + 1
             ENDIF
          ENDIF
          IF ( gradient /= 0.0_wp )  THEN
             IF ( k /= nzt )  THEN
                initial_profile(k) = initial_profile(k+1) - dzu(k+1) * gradient
             ELSE
                initial_profile(k)   = surface_value - 0.5_wp * dzu(k+1) * gradient
                initial_profile(k+1) = surface_value + 0.5_wp * dzu(k+1) * gradient
             ENDIF
          ELSE
             initial_profile(k) = initial_profile(k+1)
          ENDIF
!
!--       Avoid negative values of scalars
          IF ( initial_profile(k) < 0.0_wp )  THEN
             initial_profile(k) = 0.0_wp
          ENDIF
       ENDDO

    ENDIF

!
!-- In case of no given gradients, choose zero gradient conditions
    IF ( vertical_gradient_level(1) == -999999.9_wp )  THEN
       vertical_gradient_level(1) = 0.0_wp
    ENDIF
!
!-- Store gradient at the top boundary for possible Neumann boundary condition
    bc_top_gradient  = ( initial_profile(nzt+1) - initial_profile(nzt) ) / dzu(nzt+1)

 END SUBROUTINE init_vertical_profiles