source: LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90 @ 4545

module atke_exchange_coeff_mod

implicit none

contains 

subroutine atke_compute_km_kh(ngrid,nlay, &
                        wind_u,wind_v,temp,play,pinterf, &
                        tke,Km_out,Kh_out)

!========================================================================
! Routine that computes turbulent Km / Kh coefficients with a
! 1.5 order closure scheme (TKE) with or without stationarity assumption
!
! This parameterization has been constructed in the framework of a
! collective and collaborative workshop, 
! the so-called 'Atelier TKE (ATKE)' with
! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, T Dubos, 
! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon
!
! Main assumptions of the model :
! (1) dry atmosphere
! (2) horizontal homogeneity (Dx=Dy=0.)
!=======================================================================



USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd
USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd
USE atke_turbulence_ini_mod, ONLY : viscom, viscoh

implicit none


! Declarations:
!=============

INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
INTEGER, INTENT(IN) :: nlay ! number of vertical index  

REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_u   ! zonal velocity (m/s)
REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)


REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers

REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers
REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers

! Local variables
REAL, DIMENSION(ngrid,nlay+1) :: Km          ! Exchange coefficient for momentum at interface between layers
REAL, DIMENSION(ngrid,nlay+1) :: Kh          ! Exchange coefficient for heat flux at interface between layers
REAL, DIMENSION(ngrid,nlay)   :: theta       ! Potential temperature
REAL, DIMENSION(ngrid,nlay+1) :: l_exchange  ! Length of exchange (at interface)
REAL, DIMENSION(ngrid,nlay+1) :: z_interf    ! Altitude at the interface
REAL, DIMENSION(ngrid,nlay)   :: z_lay       ! Altitude of layers
REAL, DIMENSION(ngrid,nlay)   :: dz_interf   ! distance between two consecutive interfaces
REAL, DIMENSION(ngrid,nlay)   :: dz_lay      ! distance between two layer middles (NB: first and last are half layers)
REAL, DIMENSION(ngrid,nlay+1) :: Ri          ! Richardson's number (at interface)
REAL, DIMENSION(ngrid,nlay+1) :: Prandtl     ! Turbulent Prandtl's number (at interface)
REAL, DIMENSION(ngrid,nlay+1) :: Sm          ! Stability function for momentum (at interface)
REAL, DIMENSION(ngrid,nlay+1) :: Sh          ! Stability function for heat (at interface)

INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
REAL    :: cn,Ri0,Ri1    ! parameter for Sm stability function and Prandlt
REAL    :: preff      ! reference pressure for potential temperature calculations
REAL    :: thetam     ! mean potential temperature at interface


! Initializations:
!================

DO igrid=1,ngrid
    dz_interf(igrid,1) = 0.0
    z_interf(igrid,1) = 0.0
END DO

! Calculation of potential temperature: (if vapor -> todo virtual potential temperature)
!=====================================

preff=100000.
! The result should not depend on the choice of preff
DO ilay=1,nlay
     DO igrid = 1, ngrid
        theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
     END DO
END DO



! Calculation of altitude of layers' middle and bottom interfaces:
!=================================================================

DO ilay=2,nlay+1
    DO igrid=1,ngrid
        dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
    ENDDO
ENDDO

DO ilay=1,nlay
   DO igrid=1,ngrid 
      z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
   ENDDO
ENDDO


! Computing the mixing length:
! so far, we have neglected the effect of local stratification
!==============================================================


DO ilay=2,nlay+1
    DO igrid=1,ngrid
        l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
    ENDDO
ENDDO

! Computes the gradient Richardson's number and stability functions:
!===================================================================

! calculation of cn = Sm value at Ri=0
! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
cn=(1./sqrt(cepsilon))**(2/3)
! calculation of Ri0 such that continuity in slope of Sm at Ri=0
Ri0=2./rpi*(cinf - cn)*ric/cn
! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope


DO ilay=2,nlay
    DO igrid=1,ngrid
        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
        Ri(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) / & 
        MAX(((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & 
        ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2,1E-10) 
        
        IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
            Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn
            Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut 
        ELSE ! stable cases
            Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric))
            Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope
            IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN
               call abort_physic("atke_compute_km_kh", &
               'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)
            ENDIF
        END IF
        
        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)

    ENDDO
ENDDO

! Computing the TKE:
!===================
IF (iflag_atke == 0) THEN

! stationary solution neglecting the vertical transport of TKE by turbulence
    DO ilay=2,nlay
        DO igrid=1,ngrid
            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
            (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & 
            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) * &
            (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
        ENDDO
    ENDDO

ELSE ! TO DO
    
   call abort_physic("atke_compute_km_kh", &
        'traitement non-stationnaire de la tke pas encore prevu', 1)

END IF


! Computing eddy diffusivity coefficients:
!========================================
DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
    DO igrid=1,ngrid
        ! we add the molecular viscosity to Km,h 
        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
    END DO
END DO

! for output:
!===========
Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)

end subroutine atke_compute_km_kh


end module atke_exchange_coeff_mod