!! Fortran version of different diagnostics
! L. Fita. LMD May 2016
! gfortran module_generic.o module_ForDiagnosticsVars.o -c module_ForDiagnostics.F90
!
! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90

MODULE module_ForDiagnostics

  USE module_definitions
  USE module_generic
  USE module_ForDiagnosticsVars

  CONTAINS

!!!!!!! Calculations
! compute_cape_afwa4D: Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute 
!   CAPE, CIN, ZLFC, PLFC, LI
! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure being 
!   3rd dimension the z-dim 
! compute_cllmh3D3: Computation of low, medium and high cloudiness from a 3D CLDFRA and pressure being 
!   3rd dimension the z-dim 
! compute_cllmh: Computation of low, medium and high cloudiness
! compute_clt4D3: Computation of total cloudiness from a 4D CLDFRA being 3rd dimension the z-dim
! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim
! compute_clt: Computation of total cloudiness
! compute_psl_ecmwf: Compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
! compute_massvertint1D: Subroutine to vertically integrate a 1D variable in eta vertical coordinates
! compute_vertint1D: Subroutine to vertically integrate a 1D variable in any vertical coordinates
! compute_zint4D: Subroutine to vertically integrate a 4D variable in any vertical coordinates
! compute_zmla_generic4D: Subroutine to compute pbl-height following a generic method
! compute_zwind4D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
! compute_zwind_log4D: Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
! compute_zwindMCO3D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodolog

!!!
! Calculations
!!!

  SUBROUTINE compute_cllmh4D2(cldfra4D, pres4D, cllmh4D2, d1, d2, d3, d4)
! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA and pressure 
!   where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> cllmh(3,d1,d3,d4) 1: low, 2: medium, 3: high
! It should be properly done via an 'INTERFACE', but...

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D, pres4D
    REAL(r_k), DIMENSION(3,d1,d3,d4), INTENT(out)        :: cllmh4D2

! Local
    INTEGER                                              :: i,j,k, zdim, Ndim

!!!!!!! Variables
! cldfra4D: 4D cloud fraction values [1]
! pres4D: 4D pressure values [Pa]
! Ndim: number of dimensions of the input data
! d[1-4]: dimensions of 'cldfra4D'
! zdim: number of the vertical-dimension within the matrix
! cltlmh4D2: low, medium, high cloudiness for the 4D cldfra and d2 being 'zdim'

    fname = 'compute_cllmh4D2'
    zdim = 2
    Ndim = 4

    DO i=1, d1
      DO j=1, d3
        DO k=1, d4
          cllmh4D2(:,i,j,k) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
        END DO
      END DO
    END DO
    
    RETURN 

  END SUBROUTINE compute_cllmh4D2

  SUBROUTINE compute_cllmh3D1(cldfra3D, pres3D, cllmh3D1, d1, d2, d3)
! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA and pressure 
!   where zdim is the 1st dimension (thus, cldfra3D(d1,d2,d3) --> cllmh(3,d2,d3) 1: low, 2: medium, 3: high
! It should be properly done via an 'INTERFACE', but...

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3
    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D, pres3D
    REAL(r_k), DIMENSION(3,d2,d3), INTENT(out)           :: cllmh3D1

! Local
    INTEGER                                              :: i,j,k, zdim, Ndim

!!!!!!! Variables
! cldfra3D: 3D cloud fraction values [1]
! pres3D: 3D pressure values [Pa]
! Ndim: number of dimensions of the input data
! d[1-3]: dimensions of 'cldfra3D'
! zdim: number of the vertical-dimension within the matrix
! cltlmh3D1: low, medium, high cloudiness for the 3D cldfra and d1 being 'zdim'

    fname = 'compute_cllmh3D1'
    zdim = 1
    Ndim = 3

    DO i=1, d1
      DO j=1, d2
        cllmh3D1(:,i,j) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
      END DO
    END DO
    
    RETURN 

  END SUBROUTINE compute_cllmh3D1

  SUBROUTINE compute_cllmh(cldfra1D, cldfra2D, cldfra3D, cldfra4D, pres1D, pres2D, pres3D, pres4D,    &
    Ndim, zdim, cllmh1D, cllmh2D1, cllmh2D2, cllmh3D1, cllmh3D2, cllmh3D3, cllmh4D1, cllmh4D2,        &
    cllmh4D3, cllmh4D4, d1, d2, d3, d4)
! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D, pres1D
    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in)    :: cldfra2D, pres2D
    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D, pres3D
    REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL,                                                      &
      INTENT(in)                                         :: cldfra4D, pres4D
    REAL(r_k), DIMENSION(3), OPTIONAL, INTENT(out)       :: cllmh1D
    REAL(r_k), DIMENSION(d1,3), OPTIONAL, INTENT(out)    :: cllmh2D1
    REAL(r_k), DIMENSION(d2,3), OPTIONAL, INTENT(out)    :: cllmh2D2
    REAL(r_k), DIMENSION(d2,d3,3), OPTIONAL, INTENT(out) :: cllmh3D1
    REAL(r_k), DIMENSION(d1,d3,3), OPTIONAL, INTENT(out) :: cllmh3D2
    REAL(r_k), DIMENSION(d1,d2,3), OPTIONAL, INTENT(out) :: cllmh3D3
    REAL(r_k), DIMENSION(d2,d3,d4,3), OPTIONAL,                                                       &
      INTENT(out)                                        :: cllmh4D1
    REAL(r_k), DIMENSION(d1,d3,d4,3), OPTIONAL,                                                       &
      INTENT(out)                                        :: cllmh4D2
    REAL(r_k), DIMENSION(d1,d2,d4,3), OPTIONAL,                                                       &
      INTENT(out)                                        :: cllmh4D3
    REAL(r_k), DIMENSION(d1,d2,d3,3), OPTIONAL,                                                       &
      INTENT(out)                                        :: cllmh4D4

! Local
    INTEGER                                              :: i,j,k

!!!!!!! Variables
! cldfra[1-4]D: cloud fraction values [1]
! pres[1-4]D: pressure values [Pa]
! Ndim: number of dimensions of the input data
! d[1-4]: dimensions of 'cldfra'
! zdim: number of the vertical-dimension within the matrix
! cllmh1D: low, medium and high cloudiness for the 1D cldfra
! cllmh2D1: low, medium and high cloudiness for the 2D cldfra and d1 being 'zdim'
! cllmh2D2: low, medium and high cloudiness for the 2D cldfra and d2 being 'zdim'
! cllmh3D1: low, medium and high cloudiness for the 3D cldfra and d1 being 'zdim'
! cllmh3D2: low, medium and high cloudiness for the 3D cldfra and d2 being 'zdim'
! cllmh3D3: low, medium and high cloudiness for the 3D cldfra and d3 being 'zdim'
! cllmh4D1: low, medium and high cloudiness for the 4D cldfra and d1 being 'zdim'
! cllmh4D2: low, medium and high cloudiness for the 4D cldfra and d2 being 'zdim'
! cllmh4D3: low, medium and high cloudiness for the 4D cldfra and d3 being 'zdim'
! cllmh4D4: low, medium and high cloudiness for the 4D cldfra and d4 being 'zdim'

    fname = 'compute_cllmh'

    SELECT CASE (Ndim)
      CASE (1)
        cllmh1D = var_cllmh(cldfra1D, pres1D, d1)
      CASE (2)
        IF (zdim == 1) THEN
          DO i=1, d2
            cllmh2D1(i,:) = var_cllmh(cldfra2D(:,i), pres2D(:,i), d1)
          END DO
        ELSE IF (zdim == 2) THEN
          DO i=1, d1
            cllmh2D2(i,:) = var_cllmh(cldfra2D(:,i), pres2D(i,:), d2)
          END DO
        ELSE
          PRINT *,TRIM(ErrWarnMsg('err'))
          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
          PRINT *,'    accepted values: 1,2'
          STOP
        END IF
      CASE (3)
        IF (zdim == 1) THEN
          DO i=1, d2
            DO j=1, d3
              cllmh3D1(i,j,:) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
            END DO
          END DO
        ELSE IF (zdim == 2) THEN
          DO i=1, d1
            DO j=1, d3
              cllmh3D2(i,j,:) = var_cllmh(cldfra3D(i,:,j), pres3D(i,:,j), d2)
            END DO
          END DO
        ELSE IF (zdim == 3) THEN
          DO i=1, d1
            DO j=1, d2
              cllmh3D3(i,j,:) = var_cllmh(cldfra3D(i,j,:), pres3D(i,j,:), d3)
            END DO
          END DO
        ELSE
          PRINT *,TRIM(ErrWarnMsg('err'))
          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
          PRINT *,'    accepted values: 1,2,3'
          STOP
        END IF
      CASE (4)
        IF (zdim == 1) THEN
          DO i=1, d2
            DO j=1, d3
              DO k=1, d4
                cllmh4D1(i,j,k,:) = var_cllmh(cldfra4D(:,i,j,k), pres4D(:,i,j,k), d1)
              END DO
            END DO
          END DO
        ELSE IF (zdim == 2) THEN
          DO i=1, d1
            DO j=1, d3
              DO k=1, d4
                cllmh4D2(i,j,k,:) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
              END DO
            END DO
          END DO
        ELSE IF (zdim == 3) THEN
          DO i=1, d2
            DO j=1, d3
              DO k=1, d4
                cllmh4D3(i,j,k,:) = var_cllmh(cldfra4D(i,j,:,k), pres4D(i,j,:,k), d3)
              END DO
            END DO
          END DO
        ELSE IF (zdim == 4) THEN
          DO i=1, d1
            DO j=1, d2
              DO k=1, d3
                cllmh4D4(i,j,k,:) = var_cllmh(cldfra4D(i,j,k,:), pres4D(i,j,k,:), d4)
              END DO
            END DO
          END DO
        ELSE
          PRINT *,TRIM(ErrWarnMsg('err'))
          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
          PRINT *,'    accepted values: 1,2,3,4'
          STOP
        END IF
      CASE DEFAULT
        PRINT *,TRIM(ErrWarnMsg('err'))
        PRINT *,'  ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
        STOP
      END SELECT

    RETURN 

  END SUBROUTINE compute_cllmh

  SUBROUTINE compute_clt4D2(cldfra4D, clt4D2, d1, d2, d3, d4)
! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA
!   where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> clt(d1,d3,d4)
! It should be properly done via an 'INTERFACE', but...

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D
    REAL(r_k), DIMENSION(d1,d3,d4), INTENT(out)          :: clt4D2

! Local
    INTEGER                                              :: i,j,k, zdim, Ndim

!!!!!!! Variables
! cldfra4D: 4D cloud fraction values [1]
! Ndim: number of dimensions of the input data
! d[1-4]: dimensions of 'cldfra4D'
! zdim: number of the vertical-dimension within the matrix
! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'

    fname = 'compute_clt4D2'
    zdim = 2
    Ndim = 4

    DO i=1, d1
      DO j=1, d3
        DO k=1, d4
          clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
        END DO
      END DO
    END DO
    
    RETURN 

  END SUBROUTINE compute_clt4D2

  SUBROUTINE compute_clt3D1(cldfra3D, clt3D1, d1, d2, d3)
! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA
!   where zdim is the 1st dimension (thus, cldfra4D(d1,d2,d3) --> clt(d2,d3)
! It should be properly done via an 'INTERFACE', but...

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3
    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D
    REAL(r_k), DIMENSION(d2,d3), INTENT(out)             :: clt3D1

! Local
    INTEGER                                              :: i,j,k, zdim, Ndim

!!!!!!! Variables
! cldfra3D: 3D cloud fraction values [1]
! Ndim: number of dimensions of the input data
! d[1-3]: dimensions of 'cldfra3D'
! zdim: number of the vertical-dimension within the matrix
! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'

    fname = 'compute_clt3D1'
    zdim = 1
    Ndim = 3

    DO i=1, d2
      DO j=1, d3
        clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
      END DO
    END DO
    
    RETURN 

  END SUBROUTINE compute_clt3D1

  SUBROUTINE compute_clt(cldfra1D, cldfra2D, cldfra3D, cldfra4D, Ndim, zdim, clt1D, clt2D1, clt2D2,   &
    clt3D1, clt3D2, clt3D3, clt4D1, clt4D2, clt4D3, clt4D4, d1, d2, d3, d4)
! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D
    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in)    :: cldfra2D
    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D
    REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL,                                                      &
      INTENT(in)                                         :: cldfra4D
    REAL(r_k), OPTIONAL, INTENT(out)                     :: clt1D
    REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(out)      :: clt2D1
    REAL(r_k), DIMENSION(d2), OPTIONAL, INTENT(out)      :: clt2D2
    REAL(r_k), DIMENSION(d2,d3), OPTIONAL, INTENT(out)   :: clt3D1
    REAL(r_k), DIMENSION(d1,d3), OPTIONAL, INTENT(out)   :: clt3D2
    REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(out)   :: clt3D3
    REAL(r_k), DIMENSION(d2,d3,d4), OPTIONAL,INTENT(out) :: clt4D1
    REAL(r_k), DIMENSION(d1,d3,d4), OPTIONAL,INTENT(out) :: clt4D2
    REAL(r_k), DIMENSION(d1,d2,d4), OPTIONAL,INTENT(out) :: clt4D3
    REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL,INTENT(out) :: clt4D4

! Local
    INTEGER                                              :: i,j,k

!!!!!!! Variables
! cldfra[1-4]D: cloud fraction values [1]
! Ndim: number of dimensions of the input data
! d[1-4]: dimensions of 'cldfra'
! zdim: number of the vertical-dimension within the matrix
! clt1D: total cloudiness for the 1D cldfra
! clt2D1: total cloudiness for the 2D cldfra and d1 being 'zdim'
! clt2D2: total cloudiness for the 2D cldfra and d2 being 'zdim'
! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
! clt3D2: total cloudiness for the 3D cldfra and d2 being 'zdim'
! clt3D3: total cloudiness for the 3D cldfra and d3 being 'zdim'
! clt4D1: total cloudiness for the 4D cldfra and d1 being 'zdim'
! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
! clt4D3: total cloudiness for the 4D cldfra and d3 being 'zdim'
! clt4D4: total cloudiness for the 4D cldfra and d4 being 'zdim'

    fname = 'compute_clt'

    SELECT CASE (Ndim)
      CASE (1)
        clt1D = var_clt(cldfra1D, d1)
      CASE (2)
        IF (zdim == 1) THEN
          DO i=1, d2
            clt2D1(i) = var_clt(cldfra2D(:,i), d1)
          END DO
        ELSE IF (zdim == 2) THEN
          DO i=1, d1
            clt2D2(i) = var_clt(cldfra2D(:,i), d2)
          END DO
        ELSE
          PRINT *,TRIM(ErrWarnMsg('err'))
          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
          PRINT *,'    accepted values: 1,2'
          STOP
        END IF
      CASE (3)
        IF (zdim == 1) THEN
          DO i=1, d2
            DO j=1, d3
              clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
            END DO
          END DO
        ELSE IF (zdim == 2) THEN
          DO i=1, d1
            DO j=1, d3
              clt3D2(i,j) = var_clt(cldfra3D(i,:,j), d2)
            END DO
          END DO
        ELSE IF (zdim == 3) THEN
          DO i=1, d1
            DO j=1, d2
              clt3D3(i,j) = var_clt(cldfra3D(i,j,:), d3)
            END DO
          END DO
        ELSE
          PRINT *,TRIM(ErrWarnMsg('err'))
          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
          PRINT *,'    accepted values: 1,2,3'
          STOP
        END IF
      CASE (4)
        IF (zdim == 1) THEN
          DO i=1, d2
            DO j=1, d3
              DO k=1, d4
                clt4D1(i,j,k) = var_clt(cldfra4D(:,i,j,k), d1)
              END DO
            END DO
          END DO
        ELSE IF (zdim == 2) THEN
          DO i=1, d1
            DO j=1, d3
              DO k=1, d4
                clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
              END DO
            END DO
          END DO
        ELSE IF (zdim == 3) THEN
          DO i=1, d2
            DO j=1, d3
              DO k=1, d4
                clt4D3(i,j,k) = var_clt(cldfra4D(i,j,:,k), d3)
              END DO
            END DO
          END DO
        ELSE IF (zdim == 4) THEN
          DO i=1, d1
            DO j=1, d2
              DO k=1, d3
                clt4D4(i,j,k) = var_clt(cldfra4D(i,j,k,:), d4)
              END DO
            END DO
          END DO
        ELSE
          PRINT *,TRIM(ErrWarnMsg('err'))
          PRINT *,'  ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
          PRINT *,'    accepted values: 1,2,3,4'
          STOP
        END IF
      CASE DEFAULT
        PRINT *,TRIM(ErrWarnMsg('err'))
        PRINT *,'  ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
        STOP
      END SELECT

    RETURN 

  END SUBROUTINE compute_clt

  SUBROUTINE compute_massvertint1D(var, mutot, dz, deta, integral)
    ! Subroutine to vertically integrate a 1D variable in eta vertical coordinates

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dz
    REAL(r_k), INTENT(in)                                :: mutot
    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: var, deta
    REAL(r_k), INTENT(out)                               :: integral

! Local
    INTEGER                                              :: k

!!!!!!! Variables
! var: vertical variable to integrate (assuming kgkg-1)
! mutot: total dry-air mass in column
! dz: vertical dimension
! deta: eta-levels difference between full eta-layers

    fname = 'compute_massvertint1D'

!    integral=0.
!    DO k=1,dz
!      integral = integral + var(k)*deta(k)
!    END DO
     integral = SUM(var*deta)

    integral=integral*mutot/g

    RETURN

  END SUBROUTINE compute_massvertint1D

  SUBROUTINE compute_zint4D(var4D, dlev, zweight, d1, d2, d3, d4, int3D)
    ! Subroutine to vertically integrate a 4D variable in any vertical coordinates

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1,d2,d3,d4
    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: var4D, dlev, zweight
    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: int3D

! Local
    INTEGER                                              :: i,j,l

!!!!!!! Variables
! var4D: vertical variable to integrate
! dlev: height of layers
! zweight: weight for each level to be applied (=1. for no effect)

    fname = 'compute_zint4D'

    DO i=1,d1
      DO j=1,d2
        DO l=1,d4
          CALL compute_vertint1D(var4D(i,j,:,l),d3, dlev(i,j,:,l), zweight(i,j,:,l), &
            int3D(i,j,l))
        END DO
      END DO
    END DO

    RETURN

  END SUBROUTINE compute_zint4D

  SUBROUTINE compute_vertint1D(var, dz, deta, zweight, integral)
    ! Subroutine to vertically integrate a 1D variable in any vertical coordinates

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: dz
    REAL(r_k), DIMENSION(dz), INTENT(in)                 :: var, deta, zweight
    REAL(r_k), INTENT(out)                               :: integral

! Local
    INTEGER                                              :: k

!!!!!!! Variables
! var: vertical variable to integrate
! dz: vertical dimension
! deta: eta-levels difference between layers
! zweight: weight for each level to be applied (=1. for no effect)

    fname = 'compute_vertint1D'

!    integral=0.
!    DO k=1,dz
!      integral = integral + var(k)*deta(k)
!    END DO
    integral = SUM(var*deta*zweight)

    RETURN

  END SUBROUTINE compute_vertint1D

  SUBROUTINE compute_cape_afwa4D(ta, hur, press, zg, hgt, cape, cin, zlfc, plfc, li, parcelmethod,    &
    d1, d2, d3, d4)
! Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute CAPE, CIN, ZLFC, PLFC, LI

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4, parcelmethod
    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ta, hur, press, zg
    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: cape, cin, zlfc, plfc, li
 
! Local
    INTEGER                                              :: i, j, it
    INTEGER                                              :: ofunc

!!!!!!! Variables
! ta: air temperature [K]
! hur: relative humidity [%]
! press: air pressure [Pa]
! zg: geopotential height [gpm]
! hgt: topographical height [m]
! cape: Convective available potential energy [Jkg-1]
! cin: Convective inhibition [Jkg-1]
! zlfc: height at the Level of free convection [m]
! plfc: pressure at the Level of free convection [Pa]
! li: lifted index [1]
! parcelmethod:
!   Most Unstable = 1 (default)
!   Mean layer = 2
!   Surface based = 3

    fname = 'compute_cape_afwa4D'

    DO i=1, d1
      DO j=1, d2
        DO it=1, d4
          ofunc = var_cape_afwa1D(d3, ta(i,j,:,it), hur(i,j,:,it), press(i,j,:,it), zg(i,j,:,it),     &
            1, cape(i,j,it), cin(i,j,it), zlfc(i,j,it), plfc(i,j,it), li(i,j,it), parcelmethod)
          IF (zlfc(i,j,it) /= -1.) zlfc(i,j,it) = zlfc(i,j,it) - hgt(i,j)
        END DO
      END DO
    END DO

    RETURN

  END SUBROUTINE compute_cape_afwa4D

  SUBROUTINE compute_psl_ecmwf(ps, hgt, T, press, unpress, psl, d1, d2, d4)
! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d4
    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in)           :: ps, T, press, unpress
    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: psl
 
! Local
    INTEGER                                              :: i, j, it

!!!!!!! Variables
! ps: surface pressure [Pa]
! hgt: terrain height [m]
! T: temperature at first half-mass level [K]
! press: pressure at first full levels [Pa]
! unpress: pressure at first mass (half) levels [Pa]
! psl: sea-level pressure [Pa]

    fname = 'compute_psl_ecmwf'

    DO i=1, d1
      DO j=1, d2
        DO it=1, d4
          CALL var_psl_ecmwf(ps(i,j,it), hgt(i,j), T(i,j,it), unpress(i,j,it), press(i,j,it),         &
            psl(i,j,it))
        END DO
      END DO
    END DO

    RETURN

  END SUBROUTINE compute_psl_ecmwf

  SUBROUTINE compute_zmla_generic4D(tpot, qratio, z, hgt, zmla3D, d1, d2, d3, d4)
! Subroutine to compute pbl-height following a generic method
!    from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
!    applied also in Garcia-Diez et al., 2013, QJRMS
!   where 
!     "The technique identifies the ML height as a threshold increase of potential temperature from 
!       its minimum value within the boundary layer."
!   here applied similarly to Garcia-Diez et al. where 
!      zmla = "...first level where potential temperature exceeds the minimum potential temperature
!        reached in the mixed layer by more than 1.5 K"

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: tpot, qratio, z
    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: hgt
    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: zmla3D
 
! Local
    INTEGER                                              :: i, j, it

!!!!!!! Variables
! tpot: potential air temperature [K]
! qratio: water vapour mixing ratio [kgkg-1]
! z: height above sea level [m]
! hgt: terrain height [m]
! zmla3D: boundary layer height from surface [m]

    fname = 'compute_zmla_generic4D'

    DO i=1, d1
      DO j=1, d2
        DO it=1, d4
          CALL var_zmla_generic(d3, qratio(i,j,:,it), tpot(i,j,:,it), z(i,j,:,it), hgt(i,j),          &
            zmla3D(i,j,it))
        END DO
      END DO
    END DO

    RETURN

  END SUBROUTINE compute_zmla_generic4D

  SUBROUTINE compute_zwind4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ua, va, z
    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in)           :: uas, vas
    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: sina, cosa
    REAL(r_k), INTENT(in)                                :: zextrap
    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: uaz, vaz
 
! Local
    INTEGER                                              :: i, j, it

!!!!!!! Variables
! tpot: potential air temperature [K]
! qratio: water vapour mixing ratio [kgkg-1]
! z: height above surface [m]
! sina, cosa: local sine and cosine of map rotation [1.]
! zmla3D: boundary layer height from surface [m]

    fname = 'compute_zwind4D'

    DO i=1, d1
      DO j=1, d2
        DO it=1, d4
          CALL var_zwind(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it),       &
            sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
        END DO
      END DO
    END DO

    RETURN

  END SUBROUTINE compute_zwind4D

  SUBROUTINE compute_zwind_log4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
! Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: ua, va, z
    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in)           :: uas, vas
    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: sina, cosa
    REAL(r_k), INTENT(in)                                :: zextrap
    REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out)          :: uaz, vaz
 
! Local
    INTEGER                                              :: i, j, it

!!!!!!! Variables
! tpot: potential air temperature [K]
! qratio: water vapour mixing ratio [kgkg-1]
! z: height above surface [m]
! sina, cosa: local sine and cosine of map rotation [1.]
! zmla3D: boundary layer height from surface [m]

    fname = 'compute_zwind_log4D'

    DO i=1, d1
      DO j=1, d2
        DO it=1, d4
          CALL var_zwind_log(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it),   &
            sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
        END DO
      END DO
    END DO

    RETURN

  END SUBROUTINE compute_zwind_log4D

  SUBROUTINE compute_zwindMO3D(d1, d2, d3, ust, znt, rmol, uas, vas, sina, cosa, newz, uznew, vznew)
! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
!   NOTE: only usefull for newz < 80. m

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3
    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: ust, znt, rmol
    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: uas, vas
    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: sina, cosa
    REAL(r_k), INTENT(in)                                :: newz
    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: uznew, vznew
 
! Local
    INTEGER                                              :: i, j, it

!!!!!!! Variables
! ust: u* in similarity theory [ms-1]
! znt: thermal time-varying roughness length [m]
! rmol: Inverse of the Obukhov length [m-1]
! uas: x-component 10-m wind speed [ms-1]
! vas: y-component 10-m wind speed [ms-1]
! sina, cosa: local sine and cosine of map rotation [1.]

    fname = 'compute_zwindMO3D'

    DO i=1, d1
      DO j=1, d2
        DO it=1, d3
          CALL var_zwind_MOtheor(ust(i,j,it), znt(i,j,it), rmol(i,j,it), uas(i,j,it), vas(i,j,it),    &
            sina(i,j), cosa(i,j), newz, uznew(i,j,it), vznew(i,j,it))
        END DO
      END DO
    END DO

    RETURN

  END SUBROUTINE compute_zwindMO3D

  SUBROUTINE compute_potevap_orPM3D(d1, d2, d3, rho1, ust, uas, vas, tas, ps, qv1, potevap)
! Subroutine to compute potential evapotranspiration Penman-Monteith formulation implemented in 
!   ORCHIDEE in src_sechiba/enerbil.f90

    IMPLICIT NONE

    INTEGER, INTENT(in)                                  :: d1, d2, d3
    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: rho1, ust, uas, vas, tas, ps, qv1
    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: potevap
 
! Local
    INTEGER                                              :: i, j, it

!!!!!!! Variables
! rho1: atsmophere density at the first layer [kgm-3]
! ust: u* in similarity theory [ms-1]
! uas: x-component 10-m wind speed [ms-1]
! vas: y-component 10-m wind speed [ms-1]
! tas: 2-m atmosphere temperature [K]
! ps: surface pressure [Pa]
! qv1: 1st layer atmospheric mixing ratio [kgkg-1]
! potevap: potential evapo transpiration [kgm-2s-1]

    fname = 'compute_potevap_orPM3D'

    DO i=1, d1
      DO j=1, d2
        DO it=1, d3
          CALL var_potevap_orPM(rho1(i,j,it), ust(i,j,it), uas(i,j,it), vas(i,j,it), tas(i,j,it),     &
            ps(i,j,it), qv1(i,j,it), potevap(i,j,it))
        END DO
      END DO
    END DO

    RETURN

  END SUBROUTINE compute_potevap_orPM3D

END MODULE module_ForDiagnostics