module disvert_icosa_lmdz_mod

implicit none

contains

  SUBROUTINE disvert_icosa_lmdz(ap,bp,presnivs,presinter)
  ! Routine to compute ap(),bp() and presnivs() in the same way as done
  ! in LMDZ for planets (see LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F)

  USE prec, ONLY : rstd
  USE getin_mod, ONLY : getin
  USE grid_param, ONLY: llm
  USE earth_const, ONLY: pa,preff,scale_height
  USE abort_mod, ONLY: dynamico_abort
  USE mpipara, ONLY: is_mpi_root
  USE omp_para, ONLY: omp_in_parallel
  USE transfert_omp_mod, ONLY: bcast_omp
  USE free_unit_mod, ONLY : free_unit
  IMPLICIT NONE

  REAL(rstd),INTENT(OUT) :: ap(:)
  REAL(rstd),INTENT(OUT) :: bp(:)
  REAL(rstd),INTENT(OUT) :: presnivs(:)
  REAL(rstd),INTENT(OUT) :: presinter(:)

  INTEGER :: unit
  CHARACTER(len=255) :: filename
  INTEGER :: l,ok
  LOGICAL :: hybrid
  REAL(rstd) :: scaleheight ! atmospheric scale height (km) from file z2sig.def
  REAL(rstd) :: zsig(llm)
  REAL(rstd) :: sig(llm+1)
  REAL(rstd) :: newsig
  
    filename="z2sig.def" !file to read
    
    ! User may choose between hybrid (default) or sigma coordinates
    hybrid=.true. ! default: use hybrid coordinates
    CALL getin('hybrid',hybrid)
    ! User must also provide "pa" if in hybrid case
    if (hybrid) then
      pa=-999 ! absurd initial value to trigger error if not set
      CALL getin('pa',pa)
      IF (pa<0) THEN
        WRITE(*,*) "disvert_icosa_lmdz error: pa=",pa
        WRITE(*,*) " set a correct value for the pa parameter"
        CALL dynamico_abort("disvert_icosa_lmdz : wrong value for pa")
      ENDIF
    endif
    
!$OMP MASTER
    unit=free_unit()
    OPEN(unit,file=filename,status="old",action="read",iostat=ok)
    IF (ok/=0) THEN
      WRITE(*,*) "disvert_icosa_lmdz error: input file ",trim(filename)," not found!"
      CALL dynamico_abort("disvert_icosa_lmdz : could not open z2sig.def file" )
    ENDIF
    ! read "scaleheight" on the first line
    READ(unit,fmt=*,iostat=ok) scaleheight
    IF (ok/=0) THEN
      WRITE(*,*) "disvert_icosa_lmdz error: failed reading scaleheight"
      CALL dynamico_abort("disvert_icosa_lmdz : could not read input file")
    ENDIF
    ! While at it check that this "scaleheight" matches the "scale_height"
    ! from the def files
    IF ((1.-(scaleheight*1000.)/scale_height)>1.e-6) THEN
      WRITE(*,*) "disvert_icosa_lmdz error: significant mismatch between"
      WRITE(*,*) " scale_height (m) from def file:",scale_height
      WRITE(*,*) " and scaleheight (km) from z2sig.def:",scaleheight
      CALL dynamico_abort("disvert_icosa_lmdz : atmospheric scale height issue")
    ENDIF
    ! read zsig() line by line, starting from the surface up to model top
    DO l=1,llm
      READ(unit,fmt=*,iostat=ok) zsig(l)
      IF (ok/=0) THEN
        WRITE(*,*) "disvert_icosa_lmdz error: failed reading zsig(l) for l=",l
        CALL dynamico_abort("disvert_icosa_lmdz : could not read input file" )
      ENDIF
    ENDDO
    CLOSE(unit)
    
    ! build sig()
    sig(1)=1.
    DO l=2,llm
      sig(l)=0.5*(exp(-zsig(l)/scaleheight)+exp(-zsig(l-1)/scaleheight))
    ENDDO
    sig(llm+1)=0 ! here we asume 0 pressure at top of atmosphere

    ! build  ap() and bp() using sig(), pa and preff
    IF (hybrid) THEN  ! use hybrid coordinates
      DO l=1,llm
        CALL sig_hybrid(sig(l),pa,preff,newsig)
        bp(l)=exp(1.-1./(newsig**2))
        ap(l)=pa*(newsig-bp(l))
      ENDDO
    ELSE ! use sigma coordinates
      ap(1:llm)=0
      bp(1:llm)=sig(1:llm)
    ENDIF
    ap(llm+1)=0
    bp(llm+1)=0
!$OMP END MASTER

    IF (omp_in_parallel()) THEN
      CALL bcast_omp(ap)
      CALL bcast_omp(bp)
    ENDIF
    
    ! build presnivs(), approximative mid-layer pressures
    DO l=1,llm
      presnivs(l)=0.5*(ap(l)+bp(l)*preff+ap(l+1)+bp(l+1)*preff)
    ENDDO
    ! build presinter(), approximative inter-layer pressures
    DO l=1,llm+1
      presinter(l)=ap(l)+bp(l)*preff
    ENDDO
    
    ! tell the world about it
    IF (is_mpi_root) THEN
!$OMP MASTER
      WRITE(*,*) "Vertical discretization set by disvert_icosa_lmdz with hybrid=",hybrid
      WRITE(*,*) "ap()=",ap
      WRITE(*,*) "bp()=",bp
      WRITE(*,*) "Approximative mid-layer pressure, assuming a surface pressure preff=",preff," Pa"
      WRITE(*,*) "and approximative mid-layer height, assuming an atmospheric scale height of ",scale_height/1000," (km)"
      DO l=1,llm
        WRITE(*,*) 'PRESNIVS(',l,')=',presnivs(l),'  Z ~ ',log(preff/presnivs(l))*scale_height/1000,       &
                   ' DZ ~ ',scale_height/1000*log((ap(l)+bp(l)*preff)/ max(ap(l+1)+bp(l+1)*preff,1.e-10))
      ENDDO
!$OMP END MASTER
    ENDIF
  
    CONTAINS
  
    SUBROUTINE sig_hybrid(sig,pa,preff,newsig)
    ! F. Forget 2002
    ! Iterative procedure to find "newsig" the sigma value for a layer
    ! based on the target values provided in z2sig.def (sig)
    ! but such that:
    ! (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
    ! iterative solving of the equation uses the fact that:
    ! where exp(1-1./x**2) become << x
    ! x      exp(1-1./x**2)/x
    ! 1           1
    ! 0.68       0.5
    ! 0.5        1.E-1
    ! 0.391      1.E-2
    ! 0.333      1.E-3
    ! 0.295      1.E-4
    ! 0.269      1.E-5
    ! 0.248      1.E-6
    ! => thus on can assume  newsig = sig*preff/pa if sig*preff/pa < 0.25
    
      IMPLICIT NONE
      REAL(rstd),INTENT(in) :: sig
      REAL(rstd),INTENT(in) :: pa
      REAL(rstd),INTENT(in) :: preff
      REAL(rstd),INTENT(out) :: newsig
      REAL(rstd) :: F,x1,x2
      INTEGER :: j
    
      newsig=sig
      x1=0
      x2=1
      IF (sig.ge.1) THEN ! handle sigma=1 case
        newsig=sig
      ELSE IF (sig*preff/pa.ge.0.25) THEN
        DO J=1,9999 ! Overkill max number of iterations
          F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig
          IF (F.gt.1) THEN
            x2=newsig
            newsig=(x1+newsig)*0.5
          ELSE
            x1=newsig
            newsig=(x2+newsig)*0.5
          ENDIF
          ! Convergence test: stop if sig is approximated to within 0.01m
          IF (abs(10.*log(F)).LT.1.E-5) EXIT
        ENDDO
      ELSE ! IF (sig*preff/pa.le.0.25) THEN
        newsig=sig*preff/pa
      ENDIF
    END SUBROUTINE sig_hybrid

  END SUBROUTINE disvert_icosa_lmdz


end module disvert_icosa_lmdz_mod