module vert_prof_dissip_icosa_lmdz_mod

implicit none

contains

  subroutine vert_prof_dissip_icosa_lmdz(llm,presnivs,preff,scale_height,zvert)
  ! Routine to compute zvert(), the multiplicative coefficient
  ! for dissipation along model levels in the same way as done
  ! in LMDZ for planets (see LMDZ.COMMON/libf/dyn3d_common/inidissip.F90)
  USE prec, ONLY : rstd
  USE getin_mod, ONLY : getin
  USE abort_mod, ONLY: dynamico_abort
  USE omp_para, ONLY: is_master

  integer,intent(in) :: llm ! number of atmospheric layers
  real(rstd),intent(in) :: presnivs(:) ! approx. pressure of model layers (Pa)
  real(rstd),intent(in) :: preff ! reference surface pressure (Pa)
  real(rstd),intent(in) :: scale_height ! atmos. scale height (m)
  real(rstd),intent(out) :: zvert(:) ! multiplicative coeff. for dissipation
  
  real(rstd) :: fac_mid,fac_up
  real(rstd) :: startalt,dissip_deltaz
  real(rstd) :: dissip_hdelta,dissip_pupstart
  real(rstd) :: zz,middle,pup
  character(len=80) :: vert_prof_dissip_setup
  integer ::l
  
  ! get parameters
  vert_prof_dissip_setup="none"
  call getin("vert_prof_dissip_setup",vert_prof_dissip_setup)
  
  select case (trim(vert_prof_dissip_setup))
    case("none") ! no vertical variation
      zvert(1:llm)=1.
    
    case("altitude") ! the "Mars"/"generic" model setup
      fac_mid=3. ! multiplicative coefficient in the middle atmosphere
      call getin("dissip_fac_mid",fac_mid)
      fac_up=30. ! multiplicative coefficient in the upper atmosphere
      call getin("dissip_fac_up",fac_up)
      startalt=70. ! pseudo-altitude (km) for the middle to upper transition
      call getin("dissip_startalt",startalt)
      dissip_deltaz=30. ! size (km) of the middle-upper transition region
      call getin("dissip_deltaz",dissip_deltaz)
      
      middle=startalt+dissip_deltaz/2
      do l=1,llm
        zz= 1. - preff/presnivs(l)
        zvert(l)=fac_mid-(fac_mid-1)/(1+zz*zz)
        zvert(l)=zvert(l)*(1.0+((fac_up/fac_mid-1)*    &
                 (1-(0.5*(1+tanh(-6./                  &
                  dissip_deltaz*(scale_height/1000.*   &
                  (-log(presnivs(l)/preff))-middle)))) &
                 )))
      enddo

    case("pressure") ! the "Venus" setup
      fac_mid=2. ! multiplicative coefficient in the middle atmosphere
      call getin("dissip_fac_mid",fac_mid)
      fac_up =10. ! multiplicative coefficient in the upper atmosphere
      call getin("dissip_fac_up",fac_up)
      dissip_deltaz=10. ! size (km) of the middle-upper transition zone
      call getin("dissip_deltaz",dissip_deltaz)
      dissip_hdelta=5. ! scale height (km) in the transition region
      call getin("dissip_hdelta",dissip_hdelta)
      dissip_pupstart=1.e3 ! pressure (Pa) at which the transition starts
      call getin("dissip_pupstart",dissip_pupstart)

      Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
      do l=1,llm
        zz= 1. - preff/presnivs(l)
        zvert(l)=fac_mid-(fac_mid-1)/(1+zz*zz)
        zvert(l)=zvert(l)*(1.0+( (fac_up/fac_mid-1)  &
                 *(1-(0.5*(1+tanh(-6./dissip_deltaz* &
                 (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
      enddo
      
    case default ! none of the above
      write(*,*) "vert_prof_dissip_icosa_lmdz: Error: ",&
                 trim(vert_prof_dissip_setup),&
                 " is not a valid option for vert_prof_dissip_setup"
      call dynamico_abort("vert_prof_dissip_icosa_lmdz :"//&
        " bad value for option vert_prof_dissip_setup")
  end select ! vert_prof_dissip_setup

  if (is_master) then
    write(*,*) "vert_prof_dissip_icosa_lmdz:"
    write(*,*) "          l       presnivs(l)         zvert(l)"
    do l=1,llm
      write(*,*)l,presnivs(l),zvert(l)
    enddo
  endif

  end subroutine vert_prof_dissip_icosa_lmdz

end module vert_prof_dissip_icosa_lmdz_mod