!
! $Id: calcaerosolstrato_rrtm.f90 5605 2025-04-07 14:05:50Z lebasn $
!
SUBROUTINE calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)

  USE clesphys_mod_h
  USE iniprint_mod_h
  USE phys_state_var_mod, ONLY: tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, tau_aero_lw_rrtm
  USE phys_local_var_mod, ONLY: mdw, tausum_aero, tausum_strat, tau_strat_443, tau_strat_550, &
       tau_strat_670, tau_strat_765, tau_strat_1020, tau_strat_10um, stratomask
  USE aero_mod
  USE dimphy
!  USE temps_mod
  USE yomcst_mod_h

!USE paramet_mod_h
IMPLICIT NONE




! Variable input
  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
  LOGICAL,INTENT(IN)                     :: debut   ! le flag de l'initialisation de la physique
  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)

! Stratospheric aerosols optical properties
  REAL, DIMENSION(klon,klev,nbands_sw_rrtm) :: tau_strat, piz_strat, cg_strat
  REAL, DIMENSION(klon,klev,nwave_sw+nwave_lw) :: tau_strat_wave
  REAL, DIMENSION(klon,klev,nbands_lw_rrtm) :: tau_lw_abs_rrtm

  INTEGER k, band, wave, i
  REAL zrho, zdz

!--calculate optical properties of the aerosol size distribution from tr_seri
  tau_strat=0.0
  piz_strat=0.0
  cg_strat=0.0
  tau_strat_wave=0.0
  tau_lw_abs_rrtm=0.0
  
!-- init tau_strat vars
  tau_strat_550(:,:) =0.0
  tau_strat_1020(:,:)=0.0
  
  CALL miecalc_aer(tau_strat, piz_strat, cg_strat, tau_strat_wave, tau_lw_abs_rrtm, paprs, debut)

!--total vertical aod at the 5 SW + 1 LW wavelengths
  DO wave=1, nwave_sw+nwave_lw
    DO k=1, klev
      tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+tau_strat_wave(:,k,wave)
    ENDDO
  ENDDO

!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
  DO band=1, nbands_sw_rrtm
    !--no stratospheric aerosol in index 1
    cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,2,band)
    piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,2,band)
    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,2,band)

    !--tropospheric and stratospheric aerosol in index 2
    cg_aero_sw_rrtm(:,:,2,band) = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
                                cg_strat(:,:,band)*piz_strat(:,:,band)*tau_strat(:,:,band) ) /                              &
                                MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                            &
                                piz_strat(:,:,band)*tau_strat(:,:,band), 1.e-15 )
    piz_aero_sw_rrtm(:,:,2,band)= ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
                                piz_strat(:,:,band)*tau_strat(:,:,band) ) /                                                 &
                                MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_strat(:,:,band), 1.e-15 )
    tau_aero_sw_rrtm(:,:,2,band)= tau_aero_sw_rrtm(:,:,2,band) + tau_strat(:,:,band)
  ENDDO

  DO band=1, nbands_lw_rrtm
    !--no stratospheric aerosols in index 1
    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,2,band)
    !--tropospheric and stratospheric aerosol in index 2
    tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + tau_lw_abs_rrtm(:,:,band) 
  ENDDO

  WHERE (tau_aero_sw_rrtm .LT. 1.e-14) piz_aero_sw_rrtm=1.0
  WHERE (tau_aero_sw_rrtm .LT. 1.e-14) tau_aero_sw_rrtm=1.e-15
  WHERE (tau_aero_lw_rrtm .LT. 1.e-14) tau_aero_lw_rrtm=1.e-15

  tausum_strat(:,:)=0.0
  DO i=1,klon
  DO k=1,klev
    IF (stratomask(i,k).GT.0.5) THEN
      tausum_strat(i,1)=tausum_strat(i,1)+tau_strat_wave(i,k,1)  !--443 nm
      tausum_strat(i,2)=tausum_strat(i,2)+tau_strat_wave(i,k,2)  !--550 nm
      tausum_strat(i,3)=tausum_strat(i,3)+tau_strat_wave(i,k,3)  !--670 nm
      tausum_strat(i,4)=tausum_strat(i,4)+tau_strat_wave(i,k,4)  !--765 nm
      tausum_strat(i,5)=tausum_strat(i,5)+tau_strat_wave(i,k,5)  !--1020 nm
      tausum_strat(i,6)=tausum_strat(i,6)+tau_strat_wave(i,k,6)  !--10 um
    ENDIF
  ENDDO
  ENDDO

  DO i=1,klon
  DO k=1,klev
    zrho=pplay(i,k)/t_seri(i,k)/RD            !air density in kg/m3
    zdz=(paprs(i,k)-paprs(i,k+1))/zrho/RG     !thickness of layer in m
        tau_strat_443(i,k)=tau_strat_wave(i,k,1)/zdz
        tau_strat_550(i,k)=tau_strat_wave(i,k,2)/zdz
        tau_strat_670(i,k)=tau_strat_wave(i,k,3)/zdz
        tau_strat_765(i,k)=tau_strat_wave(i,k,4)/zdz
        tau_strat_1020(i,k)=tau_strat_wave(i,k,5)/zdz
        tau_strat_10um(i,k)=tau_strat_wave(i,k,6)/zdz
  ENDDO
  ENDDO

END SUBROUTINE calcaerosolstrato_rrtm
