MODULE calc_radiation_mod

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: calc_radiation.f90 4 2007-02-13 11:33:16Z franke $
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.6  2004/01/30 10:17:03  raasch
! Scalar lower k index nzb replaced by 2d-array nzb_2d
!
! Revision 1.1  2000/04/13 14:42:45  schroeter
! Initial revision
!
!
! Description:
! -------------
! Calculation of the vertical divergences of the long-wave radiation-fluxes
! based on the parameterization of the cloud effective emissivity
!------------------------------------------------------------------------------!

    PRIVATE
    PUBLIC calc_radiation
    
    LOGICAL, SAVE ::  first_call = .TRUE.
    REAL, SAVE    ::  sigma = 5.67E-08

    REAL, DIMENSION(:), ALLOCATABLE, SAVE ::  lwp_ground, lwp_top, &
                                              blackbody_emission

    INTERFACE calc_radiation
       MODULE PROCEDURE calc_radiation
       MODULE PROCEDURE calc_radiation_ij
    END INTERFACE calc_radiation
 
 CONTAINS


!------------------------------------------------------------------------------!
! Call for all grid points
!------------------------------------------------------------------------------!
    SUBROUTINE calc_radiation

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE indices
       USE pegrid

       IMPLICIT NONE

       INTEGER ::  i, j, k, k_help
 
       REAL :: df_p, df_m , effective_emission_up_m, effective_emission_up_p, &
               effective_emission_down_m, effective_emission_down_p,          &
               f_up_m, f_up_p, f_down_m, f_down_p, impinging_flux_at_top,     &
               temperature


!
!--    On first call, allocate temporary arrays
       IF ( first_call )  THEN
          ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1), &
                    lwp_top(nzb:nzt+1) )
          first_call = .FALSE.
       ENDIF


       DO  i = nxl, nxr
          DO  j = nys, nyn
!
!--          Compute the liquid water path (LWP) and blackbody_emission
!--          at all vertical levels
             lwp_ground(nzb) = 0.0
             lwp_top(nzt+1)  = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)

             temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
             blackbody_emission(nzb) = sigma * temperature**4.0

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

                k_help = ( nzt+nzb+1 ) - k
                lwp_ground(k)   = lwp_ground(k-1) + rho_surface * ql(k,j,i) * &
                                  dzw(k)

                lwp_top(k_help) = lwp_top(k_help+1) + &
                                  rho_surface * ql(k_help,j,i) * dzw(k_help)

                temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
                blackbody_emission(k) = sigma * temperature**4.0

             ENDDO

             lwp_ground(nzt+1) = lwp_ground(nzt) + &
                                 rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
             lwp_top(nzb)      = lwp_top(nzb+1)

             temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp * &
                                 ql(nzt+1,j,i)
             blackbody_emission(nzt+1) = sigma * temperature**4.0

!
!--          See Chlond '92, this is just a first guess
             impinging_flux_at_top = blackbody_emission(nzb) - 100.0

             DO  k = nzb_2d(j,i)+1, nzt
!
!--             Save some computational time, but this may cause load
!--             imbalances if ql is not distributed uniformly
                IF ( ql(k,j,i) /= 0.0 ) THEN
!
!--                Compute effective emissivities
                   effective_emission_up_p   = 1.0 - &
                                               EXP( -130.0 * lwp_ground(k+1) )
                   effective_emission_up_m   = 1.0 - &
                                               EXP( -130.0 * lwp_ground(k-1) )
                   effective_emission_down_p = 1.0 - &
                                               EXP( -158.0 * lwp_top(k+1) )
                   effective_emission_down_m = 1.0 - &
                                               EXP( -158.0 * lwp_top(k-1) )  

!
!--                Compute vertical long wave radiation fluxes
                   f_up_p = blackbody_emission(nzb) + &
                            effective_emission_up_p * &
                           ( blackbody_emission(k) - blackbody_emission(nzb) )

                   f_up_m = blackbody_emission(nzb) + &
                            effective_emission_up_m * &
                           ( blackbody_emission(k-1) - blackbody_emission(nzb) )

                   f_down_p = impinging_flux_at_top + &
                              effective_emission_down_p * &
                             ( blackbody_emission(k) - impinging_flux_at_top )

                   f_down_m = impinging_flux_at_top + &
                              effective_emission_down_m * &
                             ( blackbody_emission(k-1) - impinging_flux_at_top )

!
!--                Divergence of vertical long wave radiation fluxes
                   df_p = f_up_p - f_down_p
                   df_m = f_up_m - f_down_m

!
!--                Compute tendency term         
                   tend(k,j,i) = tend(k,j,i) - &
                                ( pt_d_t(k) / ( rho_surface * cp ) * &
                                  ( df_p - df_m ) / dzw(k) )

                ENDIF

             ENDDO
          ENDDO
       ENDDO

    END SUBROUTINE calc_radiation


!------------------------------------------------------------------------------!
! Call for grid point i,j
!------------------------------------------------------------------------------!
    SUBROUTINE calc_radiation_ij( i, j )

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE indices
       USE pegrid
   
       IMPLICIT NONE

       INTEGER :: i, j, k, k_help

       REAL :: df_p, df_m , effective_emission_up_m, effective_emission_up_p, &
               effective_emission_down_m, effective_emission_down_p,          &
               f_up_m, f_up_p, f_down_m, f_down_p, impinging_flux_at_top,     &
               temperature

!
!--    On first call, allocate temporary arrays
       IF ( first_call )  THEN
          ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1), &
                    lwp_top(nzb:nzt+1) )
          first_call = .FALSE.
       ENDIF

!
!--    Compute the liquid water path (LWP) and blackbody_emission
!--    at all vertical levels
       lwp_ground(nzb) = 0.0
       lwp_top(nzt+1)  = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)

       temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
       blackbody_emission(nzb) = sigma * temperature**4.0

       DO  k = nzb_2d(j,i)+1, nzt
          k_help = ( nzt+nzb+1 ) - k
          lwp_ground(k)   = lwp_ground(k-1) + rho_surface * ql(k,j,i) * dzw(k)

          lwp_top(k_help) = lwp_top(k_help+1) + &
                            rho_surface * ql(k_help,j,i) * dzw(k_help)

          temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
          blackbody_emission(k) = sigma * temperature**4.0

       ENDDO
       lwp_ground(nzt+1) = lwp_ground(nzt) + &
                           rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
       lwp_top(nzb)      = lwp_top(nzb+1)

       temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp * &
                           ql(nzt+1,j,i)
       blackbody_emission(nzt+1) = sigma * temperature**4.0

!
!--    See Chlond '92, this is just a first guess
       impinging_flux_at_top = blackbody_emission(nzb) - 100.0

       DO  k = nzb_2d(j,i)+1, nzt
!
!--       Store some computational time,
!--       this may cause load imbalances if ql is not distributed uniformly
          IF ( ql(k,j,i) /= 0.0 ) THEN
!
!--          Compute effective emissivities
             effective_emission_up_p   = 1.0 - &
                                         EXP( -130.0 * lwp_ground(k+1) )
             effective_emission_up_m   = 1.0 - &
                                         EXP( -130.0 * lwp_ground(k-1) )
             effective_emission_down_p = 1.0 - &
                                         EXP( -158.0 * lwp_top(k+1) )
             effective_emission_down_m = 1.0 - &
                                         EXP( -158.0 * lwp_top(k-1) )  
             
!
!--          Compute vertical long wave radiation fluxes
             f_up_p = blackbody_emission(nzb) + effective_emission_up_p * &
                     ( blackbody_emission(k) - blackbody_emission(nzb) )

             f_up_m = blackbody_emission(nzb) + effective_emission_up_m * &
                     ( blackbody_emission(k-1) - blackbody_emission(nzb) )

             f_down_p = impinging_flux_at_top + effective_emission_down_p * &
                       ( blackbody_emission(k) - impinging_flux_at_top )

             f_down_m = impinging_flux_at_top + effective_emission_down_m * &
                       ( blackbody_emission(k-1) - impinging_flux_at_top )

!
!-           Divergence of vertical long wave radiation fluxes
             df_p = f_up_p - f_down_p
             df_m = f_up_m - f_down_m

!
!--          Compute tendency term         
             tend(k,j,i) = tend(k,j,i) - ( pt_d_t(k) / ( rho_surface * cp ) * &
                                         ( df_p - df_m ) / dzw(k) )

          ENDIF

       ENDDO

    END SUBROUTINE calc_radiation_ij

 END MODULE calc_radiation_mod