MODULE nudge_mod

!--------------------------------------------------------------------------------!
! This file is part of PALM.
!
! 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 1997-2012  Leibniz University Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
!
! Former revisions:
! -----------------
! $Id: nudging.f90 1269 2013-12-12 09:49:00Z fricke $
!
! 1268 2013-12-12 09:47:53Z heinze
! bugfix: argument of calc_mean_profile corrected
!
! 1251 2013-11-07 08:14:30Z heinze
! bugfix: calculate dtm and dtp also in vector version
!
! 1249 2013-11-06 10:45:47Z heinze
! remove call of user module
! reformatting
!
! 1241 2013-10-30 11:36:58Z heinze
! Initial revision
!
! Description:
! ------------
! Nudges u, v, pt and q to given profiles on a relaxation timescale tnudge. 
! Profiles are read in from NUDGIN_DATA. Code is based on Neggers et al. (2012) 
! and also part of DALES and UCLA-LES.
!--------------------------------------------------------------------------------!

    PRIVATE
    PUBLIC init_nudge, nudge 
    SAVE

    INTERFACE nudge
       MODULE PROCEDURE nudge
       MODULE PROCEDURE nudge_ij
    END INTERFACE nudge

 CONTAINS

    SUBROUTINE init_nudge

       USE arrays_3d
       USE control_parameters
       USE cpulog
       USE indices
       USE interfaces
       USE pegrid

       IMPLICIT NONE

       INTEGER :: finput = 90, ierrn, k, t 

       CHARACTER(1) :: hash
       REAL :: highheight, highqnudge, highptnudge, highunudge, highvnudge, &
               highwnudge, hightnudge
       REAL :: lowheight, lowqnudge, lowptnudge, lowunudge, lowvnudge, &
               lowwnudge, lowtnudge
       REAL :: fac

       ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
                 tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge),  &
                 vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge)  )

       ALLOCATE( timenudge(0:ntnudge) )


       ptnudge = 0.0; qnudge = 0.0; tnudge = 0.0; unudge = 0.0
       vnudge = 0.0; wnudge = 0.0; timenudge = 0.0

       t = 0
       OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
              FORM='FORMATTED', IOSTAT=ierrn )

       IF ( ierrn /= 0 )  THEN
          message_string = 'file NUDGING_DATA does not exist'
          CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
       ENDIF

       ierrn = 0

 rloop:DO
          t = t + 1
          hash = "#"
          ierrn = 1 ! not zero        
!
!--       Search for the next line consisting of "# time", 
!--       from there onwards the profiles will be read
          DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) ) 
          
            READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(t)
            IF ( ierrn < 0 )  EXIT rloop

          ENDDO

          ierrn = 0
          READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge,   &
                                           lowvnudge, lowwnudge , lowptnudge, &
                                           lowqnudge

          IF ( ierrn /= 0 )  THEN
             message_string = 'errors in file NUDGING_DATA'
             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
          ENDIF

          ierrn = 0
          READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge,   &
                                           highvnudge, highwnudge , highptnudge, &
                                           highqnudge

          IF ( ierrn /= 0 )  THEN
             message_string = 'errors in file NUDGING_DATA'
             CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
          ENDIF

          DO  k = nzb, nzt+1
             IF ( highheight < zu(k) )  THEN
                lowheight  = highheight
                lowtnudge  = hightnudge
                lowunudge  = highunudge
                lowvnudge  = highvnudge
                lowwnudge  = highwnudge
                lowptnudge = highptnudge
                lowqnudge  = highqnudge
 
                ierrn = 0
                READ ( finput, *, IOSTAT=ierrn )  highheight , hightnudge , &
                                                  highunudge , highvnudge , &
                                                  highwnudge , highptnudge, &
                                                  highqnudge
                IF (ierrn /= 0 )  THEN
                   message_string = 'errors in file NUDGING_DATA'
                   CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
                ENDIF
             ENDIF

!
!--          Interpolation of prescribed profiles in space 

             fac = ( highheight - zu(k) ) / ( highheight - lowheight )

             tnudge(k,t)  = fac * lowtnudge + ( 1 - fac ) * hightnudge
             unudge(k,t)  = fac * lowunudge + ( 1 - fac ) * highunudge
             vnudge(k,t)  = fac * lowvnudge + ( 1 - fac ) * highvnudge
             wnudge(k,t)  = fac * lowwnudge + ( 1 - fac ) * highwnudge
             ptnudge(k,t) = fac * lowptnudge + ( 1 - fac ) * highptnudge
             qnudge(k,t)  = fac * lowqnudge + ( 1 - fac ) * highqnudge
          ENDDO

       ENDDO rloop

       CLOSE ( finput )

!
!--    Prevent nudging if nudging profiles exhibt too small values
!--    not used so far
       lptnudge  = ANY( ABS( ptnudge ) > 1e-8 )
       lqnudge   = ANY( ABS( qnudge ) > 1e-8 )
       lunudge   = ANY( ABS( unudge ) > 1e-8 )
       lvnudge   = ANY( ABS( vnudge ) > 1e-8 )
       lwnudge   = ANY( ABS( wnudge ) > 1e-8 )

    END SUBROUTINE init_nudge

!------------------------------------------------------------------------------!
! Call for all grid points
!------------------------------------------------------------------------------!
    SUBROUTINE nudge ( time, prog_var )

       USE arrays_3d
       USE buoyancy_mod
       USE control_parameters
       USE cpulog
       USE indices
       USE interfaces
       USE pegrid
       USE statistics 

       IMPLICIT NONE

       CHARACTER (LEN=*) ::  prog_var

       REAL :: currtnudge, dtm, dtp, time

       INTEGER ::  i, j, k, t


       t = 1
       DO WHILE ( time > timenudge(t) )
         t = t+1
       ENDDO
       IF ( time /= timenudge(1) ) THEN
         t = t-1
       ENDIF

       dtm = ( time - timenudge(t) ) / ( timenudge(t+1) - timenudge(t) )
       dtp = ( timenudge(t+1) - time ) / ( timenudge(t+1) - timenudge(t) )

       SELECT CASE ( prog_var )

          CASE ( 'u' )

             CALL calc_mean_profile( u, 1, 'nudging' )

             DO  i = nxl, nxr
                DO  j = nys, nyn
                   DO  k = nzb_u_inner(j,i)+1, nzt

                      currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
                                               tnudge(k,t+1) * dtm )


                      tend(k,j,i) = tend(k,j,i) - ( hom(k,1,1,0) - &
                                    ( unudge(k,t) * dtp + &
                                       unudge(k,t+1) * dtm ) ) / currtnudge
                   ENDDO
                ENDDO
            ENDDO

          CASE ( 'v' )

             CALL calc_mean_profile( v, 2, 'nudging' )

             DO  i = nxl, nxr
                DO  j = nys, nyn
                   DO  k = nzb_u_inner(j,i)+1, nzt

                      currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
                                               tnudge(k,t+1) * dtm )


                      tend(k,j,i) = tend(k,j,i) - ( hom(k,1,2,0) - &
                                    ( vnudge(k,t) * dtp + &
                                       vnudge(k,t+1) * dtm ) ) / currtnudge
                   ENDDO
                ENDDO
            ENDDO

          CASE ( 'pt' )

             CALL calc_mean_profile( pt, 4, 'nudging' )

             DO  i = nxl, nxr
                DO  j = nys, nyn
                   DO  k = nzb_u_inner(j,i)+1, nzt

                      currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
                                               tnudge(k,t+1) * dtm )


                      tend(k,j,i) = tend(k,j,i) - ( hom(k,1,4,0) - &
                                    ( ptnudge(k,t) * dtp + &
                                       ptnudge(k,t+1) * dtm ) ) / currtnudge
                   ENDDO
                ENDDO
            ENDDO

          CASE ( 'q' )

             CALL calc_mean_profile( q, 41, 'nudging' )

             DO  i = nxl, nxr
                DO  j = nys, nyn
                   DO  k = nzb_u_inner(j,i)+1, nzt

                      currtnudge = MAX( dt_3d, tnudge(k,t) * dtp +  &
                                               tnudge(k,t+1) * dtm )


                      tend(k,j,i) = tend(k,j,i) - ( hom(k,1,41,0) - &
                                    ( qnudge(k,t) * dtp + &
                                       qnudge(k,t+1) * dtm ) ) / currtnudge
                   ENDDO
                ENDDO
            ENDDO

          CASE DEFAULT
             message_string = 'unknown prognostic variable "' // prog_var // '"'
             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )

       END SELECT

    END SUBROUTINE nudge


!------------------------------------------------------------------------------!
! Call for grid point i,j
!------------------------------------------------------------------------------!

    SUBROUTINE nudge_ij( i, j, time, prog_var )

       USE arrays_3d
       USE buoyancy_mod
       USE control_parameters
       USE cpulog
       USE indices
       USE interfaces
       USE pegrid
       USE statistics 

       IMPLICIT NONE

       CHARACTER (LEN=*) ::  prog_var

       REAL :: currtnudge, dtm, dtp, time

       INTEGER ::  i, j, k, t


       t = 1
       DO WHILE ( time > timenudge(t) )
         t = t+1
       ENDDO
       IF ( time /= timenudge(1) )  THEN
         t = t-1
       ENDIF

       dtm = ( time - timenudge(t) ) / ( timenudge(t+1) - timenudge(t) )
       dtp = ( timenudge(t+1) - time ) / ( timenudge(t+1) - timenudge(t) )

       SELECT CASE ( prog_var )

          CASE ( 'u' )

             CALL calc_mean_profile( u, 1, 'nudging' )

             DO  k = nzb_u_inner(j,i)+1, nzt

                currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )

                tend(k,j,i) = tend(k,j,i) - ( hom(k,1,1,0) - &
                               ( unudge(k,t) * dtp + &
                                 unudge(k,t+1) * dtm ) ) / currtnudge
             ENDDO

          CASE ( 'v' )

             CALL calc_mean_profile( v, 2, 'nudging' )

             DO  k = nzb_u_inner(j,i)+1, nzt

                currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )

                tend(k,j,i) = tend(k,j,i) - ( hom(k,1,2,0) - &
                               ( vnudge(k,t) * dtp + &
                                 vnudge(k,t+1) * dtm ) ) / currtnudge
             ENDDO


          CASE ( 'pt' )

             CALL calc_mean_profile( pt, 4, 'nudging' )

             DO  k = nzb_u_inner(j,i)+1, nzt

                currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )

                tend(k,j,i) = tend(k,j,i) - ( hom(k,1,4,0) - &
                               ( ptnudge(k,t) * dtp + &
                                 ptnudge(k,t+1) * dtm ) ) / currtnudge
             ENDDO


          CASE ( 'q' )

             CALL calc_mean_profile( q, 41, 'nudging' )

             DO  k = nzb_u_inner(j,i)+1, nzt

                currtnudge = MAX( dt_3d, tnudge(k,t) * dtp + tnudge(k,t+1) * dtm )

                tend(k,j,i) = tend(k,j,i) - ( hom(k,1,41,0) - &
                               ( qnudge(k,t) * dtp + &
                                 qnudge(k,t+1) * dtm ) ) / currtnudge
             ENDDO

          CASE DEFAULT
             message_string = 'unknown prognostic variable "' // prog_var // '"'
             CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )

       END SELECT


    END SUBROUTINE nudge_ij

 END MODULE nudge_mod