module test06_mtckd_wrapper

!!$  implicit none

  use vtype_module

  private

  public :: call_mtckd_wrapper

contains

  !----------------------------------------------------------------------------


  !----------------------------------------------------------------------------

  subroutine call_mtckd_wrapper( &
    & Arg_PAVE, Arg_TAVE, Arg_VMRH2O, Arg_xlength, &
    & MolNum, &
    & Arg_VMRCO2, Arg_VMRO3, Arg_VMRO2, Arg_VMRN2, Arg_VMROTHERS, &
    & Arg_v1abs, Arg_v2abs, Arg_dvabs, &
    & Arg_arraysize, &
    & a_absrb, &
    & Arg_NAbsrb &
    & )

    implicit real*8 (v)

    real(DP), intent(in) :: Arg_PAVE
    real(DP), intent(in) :: Arg_TAVE
    real(DP), intent(in) :: Arg_VMRH2O
    real(DP), intent(in) :: Arg_xlength
    integer , intent(in) :: MolNum
    real(DP), intent(in) :: Arg_VMRCO2
    real(DP), intent(in) :: Arg_VMRO3
    real(DP), intent(in) :: Arg_VMRO2
    real(DP), intent(in) :: Arg_VMRN2
    real(DP), intent(in) :: Arg_VMROTHERS
    real(DP), intent(in) :: Arg_v1abs
    real(DP), intent(in) :: Arg_v2abs
    real(DP), intent(in) :: Arg_dvabs
    integer , intent(in) :: Arg_arraysize
    real(DP), intent(out) :: a_absrb (Arg_arraysize)
    integer , intent(out) :: Arg_NAbsrb


    ! Local variables
    !
    logical :: Flag_Output
    integer :: Arg_ipr
    integer :: Arg_ipu
    integer :: Arg_jrad
    real(DP) :: a_radfld(Arg_arraysize)
    real(DP) :: a_csh2o (Arg_arraysize)
    real(DP) :: a_cfh2o (Arg_arraysize)
    integer  :: Arg_NH2OArray
    real(DP) :: Arg_v1h
    real(DP) :: Arg_dvh

    logical  :: Arg_Flag_h2o
    logical  :: Arg_Flag_co2
    logical  :: Arg_Flag_o3
    logical  :: Arg_Flag_o2
    logical  :: Arg_Flag_n2
    logical  :: Arg_Flag_rayl


    Arg_Flag_h2o  = .false.
    Arg_Flag_co2  = .false.
    Arg_Flag_o3   = .false.
    Arg_Flag_o2   = .false.
    Arg_Flag_n2   = .false.
    Arg_Flag_rayl = .false.

    select case ( MolNum )
    case ( -100 )
      Arg_Flag_rayl = .true.
    case ( 1 )
      Arg_Flag_h2o  = .true.
    case ( 2 )
      Arg_Flag_co2  = .true.
    case ( 3 )
      Arg_Flag_o3   = .true.
    case ( 7 )
      Arg_Flag_o2   = .true.
    case ( 22 )
      Arg_Flag_n2   = .true.
    case default
      write( 6, * ) 'Unexpected molecular number: ', MolNum
      stop
    end select

    Arg_jrad = 1


    Flag_Output = .false.

    Arg_ipr = 66
    Arg_ipu = 7


    if ( Flag_Output ) then
      open (Arg_ipr,file='CNTNM.OPTDPT')
      open (Arg_ipu,file='WATER.COEF')
    end if


    call call_mtckd( &
      & Flag_Output, &
      & Arg_ipr, Arg_ipu, &
      & Arg_jrad, &
      & Arg_PAVE, Arg_TAVE, Arg_VMRH2O, Arg_xlength, &
      & Arg_VMRCO2, Arg_VMRO3, Arg_VMRO2, Arg_VMRN2, Arg_VMROTHERS, &
      & Arg_Flag_h2o, Arg_Flag_co2, Arg_Flag_o3, Arg_Flag_o2, Arg_Flag_n2, Arg_Flag_rayl, &
      & Arg_v1abs, Arg_v2abs, Arg_dvabs, &
      & Arg_arraysize, &
      & a_absrb, &
      & Arg_NAbsrb, &
      & a_radfld, a_csh2o, a_cfh2o, &
      & Arg_NH2OArray, &
      & Arg_v1h, Arg_dvh &
      & )

    if ( Flag_Output ) then
      do i = 1, Arg_NAbsrb
        vi = Arg_v1abs + real(i-1) * Arg_dvabs
100     WRITE (Arg_ipr, 910) VI, a_ABSRB(I)
      end do
910   FORMAT(F10.3,1P,E12.3)

      WRITE (Arg_ipu,920) Arg_tave
920   FORMAT(//,' self and foreign water vapor continuum coefficients ',/,&
        'for  ',f8.2,'K - ',   //,&
        ' the self-continuum scales as ( Rself/Ro ) ',/,&
        ' the foreign continuum scales as ( (Rtot-Rself)/Ro ) ',/,&
        ' where R is the density rho [ R = (P/Po)*(To/T) ]. ',//,&
        10x,'     without radiation field:  ',&
        10x,'       with radiation field:   ',      /,&
        10x,'      self         foreign     ',&
        10x,'      self         foreign     ',      /,&
        '    cm-1  ',&
        '       1/(cm-1 molec/cm**2)    ',&
        10x,'       1/(molec/cm**2) '        ,//)

      do i=1,Arg_NH2OArray
        vi=Arg_v1h+float(i-1)*Arg_dvh
        if (vi.ge.Arg_v1abs .and. vi.le.Arg_v2abs) then
          radfld= a_radfld(i)
          csh2or= a_csh2o(i) * radfld
          cfh2or= a_cfh2o(i) * radfld
          write (Arg_ipu,930) vi, a_csh2o(i), a_cfh2o(i), csh2or, cfh2or
930       format(f10.2, 1p, 2e15.4,10x, 1p, 2e15.4)
        endif
      end do

    end if

  end subroutine call_mtckd_wrapper

  !----------------------------------------------------------------------------

  subroutine call_mtckd( &
    & Flag_Output, &
    & Arg_ipr, Arg_ipu, &
    & Arg_jrad, &
    & Arg_PAVE, Arg_TAVE, Arg_VMRH2O, Arg_xlength, &
    & Arg_VMRCO2, Arg_VMRO3, Arg_VMRO2, Arg_VMRN2, Arg_VMROTHERS, &
    & Arg_Flag_h2o, Arg_Flag_co2, Arg_Flag_o3, Arg_Flag_o2, Arg_Flag_n2, Arg_Flag_rayl, &
    & Arg_v1abs, Arg_v2abs, Arg_dvabs, &
    & Arg_arraysize, &
    & a_absrb, &
    & Arg_NAbsrb, &
    & a_radfld, a_csh2o, a_cfh2o, &
    & Arg_NH2OArray, &
    & Arg_v1h, Arg_dvh &
    & )

    implicit real*8 (v)

    logical , intent(in) :: Flag_Output
    integer, intent(in) :: Arg_ipr
    integer, intent(in) :: Arg_ipu
    integer, intent(in) :: Arg_jrad
    real(DP), intent(in) :: Arg_PAVE
    real(DP), intent(in) :: Arg_TAVE
    real(DP), intent(in) :: Arg_VMRH2O
    real(DP), intent(in) :: Arg_xlength
    real(DP), intent(in) :: Arg_VMRCO2
    real(DP), intent(in) :: Arg_VMRO3
    real(DP), intent(in) :: Arg_VMRO2
    real(DP), intent(in) :: Arg_VMRN2
    real(DP), intent(in) :: Arg_VMROTHERS
    logical , intent(in) :: Arg_Flag_h2o
    logical , intent(in) :: Arg_Flag_co2
    logical , intent(in) :: Arg_Flag_o3
    logical , intent(in) :: Arg_Flag_o2
    logical , intent(in) :: Arg_Flag_n2
    logical , intent(in) :: Arg_Flag_rayl
    real(DP), intent(in) :: Arg_v1abs
    real(DP), intent(in) :: Arg_v2abs
    real(DP), intent(in) :: Arg_dvabs
    integer , intent(in) :: Arg_arraysize
    real(DP), intent(out) :: a_absrb (Arg_arraysize)
    integer , intent(out) :: Arg_NAbsrb
    real(DP), intent(out) :: a_radfld(Arg_arraysize)
    real(DP), intent(out) :: a_csh2o (Arg_arraysize)
    real(DP), intent(out) :: a_cfh2o (Arg_arraysize)
    integer , intent(out) :: Arg_NH2OArray
    real(DP), intent(out) :: Arg_v1h
    real(DP), intent(out) :: Arg_dvh


    !
      COMMON /ABSORB/ V1ABS,V2ABS,DVABS,NPTABS,ABSRB(5050)                
!
      COMMON /CVRCNT/ HNAMCNT,HVRCNT
!
      COMMON /FILHDR/ XID(10),SECANT,PAVE,TAVE,HMOLID(60),XALTZ(4),&
                      WK(60),PZL,PZU,TZL,TZU,WBROAD,DV ,V1 ,V2 ,TBOUND,&
                      EMISIV,FSCDID(17),NMOL,LAYER ,YI1,YID(10),LSTWDF    
!
      Common /share/ HOLN2
!
      COMMON /CONSTS/ PI,PLANCK,BOLTZ,CLIGHT,AVOGAD,ALOSMT,GASCON,&
                      RADCN1,RADCN2,GRAV,CPDAIR,AIRMWT,SECDY
      COMMON /LAMCHN/ ONEPL,ONEMI,EXPMIN,ARGMIN                           
      COMMON /XCONT/  V1C,V2C,DVC,NPTC,C(6000)
!
!********************************************
      COMMON /cnth2o/ V1h,V2h,DVh,NPTh,Ch(5050),csh2o(5050),cfh2o(5050)
!********************************************
!
      COMMON /IFIL/ IRD,IPRcnt,IPU,NOPR,NFHDRF,NPHDRF,NFHDRL,NPHDRL,&
                    NLNGTH,KFILE,KPANEL,LINFIL,NFILE,IAFIL,IEXFIL,&
                    NLTEFL,LNFIL4,LNGTH4                                  

      common /cntscl/ XSELF,XFRGN,XCO2C,XO3CN,XO2CN,XN2CN,XRAYL

!------------------------------------
! for analytic derivative calculation
! note: ipts  = same dimension as ABSRB
!       ipts2 = same dimension as C
      parameter (ipts=5050,ipts2=6000)
      common /CDERIV/ icflg,iuf,v1absc,v2absc,dvabsc,nptabsc,&
          dqh2oC(ipts),dTh2oC(ipts),dUh2o,&
          dqco2C(ipts),dTco2C(ipts),&
          dqo3C(ipts),dTo3C(ipts),&
          dqo2C(ipts),dTo2C(ipts),&
          dqn2C(ipts),dTn2C(ipts)

      real ctmp(ipts2),cself(ipts),cforeign(ipts),ch2o(ipts2)
      real ctmp2(ipts2),ctmp3(ipts2),ctmp4(ipts),ctmp5(ipts)

!------------------------------------
!
      dimension xcnt(7)
!
      ! following is defined above
      !   common /cntscl/ XSELF,XFRGN,XCO2C,XO3CN,XO2CN,XN2CN,XRAYL
      ! so,
      ! xcnt(1) = xself
      ! xcnt(2) = xfrgn
      ! xcnt(3) = xco2c
      ! xcnt(4) = xo3cn
      ! xcnt(5) = xo2cn
      ! xcnt(6) = xn2cn
      ! xcnt(7) = xrayl
      !
      ! xco2c : number of particles for CO2 continuum absorption
      !       : the value is normalized to that of the present atmosphere
      ! x*cn  : the same as xco2c, but for * species
      ! xrayl : number of particles affecting Rayleigh scattering
      !       : the value is normalized to that of the present atmosphere
      !
      ! I cannot understand the difference between WK(*) and x*cn.
      !
      equivalence (xself,xcnt(1))
!
      CHARACTER*18 HNAMCNT,HVRCNT
!                                                                         F00100
      CHARACTER*8      XID,       HMOLID,      YID
      REAL*8               SECANT,       XALTZ
!
      character*8 holn2
!                                                                         F00120
      DATA XLOSMT/2.68675E+19/
!
      RADCN1=2.*PLANCK*CLIGHT*CLIGHT*1.E-07
      RADCN2=PLANCK*CLIGHT/BOLTZ
!
      icflg = -999
!
      do i=1,7
        xcnt(i)=1.
      end do
      if ( .not. Arg_Flag_h2o  ) then
        xcnt(1) = 0.0d0
        xcnt(2) = 0.0d0
      end if
      if ( .not. Arg_Flag_co2  ) xcnt(3) = 0.0d0
      if ( .not. Arg_Flag_o3   ) xcnt(4) = 0.0d0
      if ( .not. Arg_Flag_o2   ) xcnt(5) = 0.0d0
      if ( .not. Arg_Flag_n2   ) xcnt(6) = 0.0d0
      if ( .not. Arg_Flag_rayl ) xcnt(7) = 0.0d0

      do 2, i=1,5050
         absrb(i)=0.
2     continue

      do 3, i=1,60
         wk(i)=0.
 3    continue
!
      ird = 55
      ipr = Arg_ipr
      ipu = Arg_ipu
!
   if ( Flag_Output ) then
      print *
      print *, '  This version is limited to 5000 values  '
    end if
!
!   THIS PROGRAM CALCULATES THE CONTINUUM OPTICAL DEPTH
!         FOR AN HOMOGENEOUS LAYER
!
!   THE FOLLOWING QUANTITIES MUST BE SPECIFIED: 
!
!          PRESSURE                   PAVE (MB)
!
!          TEMPERATURE                TAVE ( K)
!
!          COLUMN AMOUNT
!            NITROGEN                 WN2    (MOLEC/CM**2)
!            OXYGEN                   WK(7)  (MOLEC/CM**2)
!            CARBON DIOXIDE           WK(2)  (MOLEC/CM**2)
!            WATER VAPOR              WK(1)  (MOLEC/CM**2)
!
!          NUMBER OF MOLECULES        NMOL
!
!          BEGINNING WAVENUMBER       V1ABS (CM-1)
!
!          ENDING WAVENUMBER          V2ABS (CM-1)
!
!          SAMPLING INTERVAL          DVABS (CM-1)
!
!          NUMBER OF VALUES           NPTABS
!
!
!   THE RESULTS ARE IN ARRAY ABSORB
!
!   NOTE THAT FOR AN ATMOSPHERIC LAYER: 
!
!            WTOT   = XLOSMT * (PAVE/1013) * (273/TAVE) * (PATH LENGTH)
!
!            WBROAD = the column amount for all species not explicitly provided
!
!            WK(M)  = (VOLUME MIXING RATIO) * (COLUMN OF DRY AIR)
!
!
      iprcnt = ipr
                   CALL PRCNTM
      iprcnt = ipu
                   CALL PRCNTM
!
!   THE FOLLOWING IS AN EXAMPLE FOR A ONE CM PATH (SEE CNTNM.OPTDPT FOR RESULTS)
!

      PAVE = Arg_PAVE
      TAVE = Arg_TAVE
!
      VMRH2O = Arg_VMRH2O
!
      xlength = Arg_xlength

   if ( Flag_Output ) then
      print *,&
       'Pressure (mb), Temperature (K), Path Length (cm),    VMR H2O'

      print 911, pave,tave,xlength,vmrh2o
911   format(1x,f13.6,f17.4,f18.4,f12.8)
    end if
!
!     It may be preferable to specifiy the water column directly!
!
      WTOT = XLOSMT*(PAVE/1013.)*(273./TAVE)* xlength
!
      W_dry = WTOT * (1.-VMRH2O)
!
!     ww is column of dry air;  vol. mix. ratios are based on dry air
!
! argon:
!!$      WA     = 0.009     * W_dry
!!$      WA     = Arg_VMRAr     * W_dry
      WA     = Arg_VMROTHERS * W_dry
! nitrogen:
!!$      WN2    = 0.78      * W_dry
      WN2    = Arg_VMRN2     * W_dry
! oxygen:
!!$      WK(7)  = 0.21      * W_dry
      WK(7)  = Arg_VMRO2      * W_dry

! carbon dioxide:
!!$      WK(2)  = 345.E-06  * W_dry
      WK(2)  = Arg_VMRCO2  * W_dry

! ozone:
      WK(3)  = Arg_VMRO3  * W_dry

! water vapor:
      if (abs(vmrh2o-1.) .lt. 1.e-05) then
         wk(1) = wtot
      else
         WK(1) = VMRH2O * W_dry
      endif
!
      WBROAD=WN2+WA
!
      NMOL = 7
!

!!$      V1ABS =    0.
!!$      V2ABS = 10000.
!!$
!!$      DVABS =    2.

      V1ABS = Arg_v1abs
      V2ABS = Arg_v2abs
      DVABS = Arg_dvabs
! ..........................................................
!      write (*,*) '  v1abs,  v2abs,  dvabs  '
!      read  (*,*)    v1abs,  v2abs,  dvabs
! ..........................................................

      NPTABS =  1. + (v2abs-v1abs)/dvabs

      do 85 i=1,nptabs
         absrb(i) =0.
 85   continue

   if ( Flag_Output ) then
!c
      WRITE (IPR,970) PAVE,TAVE
      WRITE (IPR,975) (HMOLID(I),I=1,7),HOLN2
      WRITE (IPR,980) (WK(M),M=1,7),WBROAD
!
      WRITE (IPu,970) PAVE,TAVE
      WRITE (IPu,975) (HMOLID(I),I=1,7),HOLN2
      WRITE (IPu,980) (WK(M),M=1,7),WBROAD
!
  970 FORMAT (/,29x, 'P(MB)',7X,'T(K)', //,23x,0P,F12.3,F9.2)
  975 FORMAT (/, 9X,'MOLECULAR AMOUNTS (MOL/CM**2) BY LAYER ',//,&
               8(1X,A6,3X))
980   FORMAT (/,1P,8E10.3,//)

    end if
!
      jrad=1
      jrad = Arg_jrad
!
      v1 = v1abs
      v2 = v2abs
!
      CALL CONTNM(JRAD)
!
      Arg_NAbsrb = nptabs
      do i = 1, nptabs
        a_absrb(i) = absrb(i)
      end do

!!$      DO 100 I=1,NPTABS
!!$        VI=V1ABS+FLOAT(I-1)*DVABS
!!$100   WRITE (ipr, 910) VI, ABSRB(I)
!!$910   FORMAT(F10.3,1P,E12.3)
!
!!$      WRITE (7,920) tave
!!$  920 FORMAT(//,' self and foreign water vapor continuum coefficients ',/,&
!!$             'for  ',f8.2,'K - ',   //,&
!!$             ' the self-continuum scales as ( Rself/Ro ) ',/,&
!!$             ' the foreign continuum scales as ( (Rtot-Rself)/Ro ) ',/,&
!!$             ' where R is the density rho [ R = (P/Po)*(To/T) ]. ',//,&
!!$         10x,'     without radiation field:  ',&
!!$         10x,'       with radiation field:   ',      /,&
!!$         10x,'      self         foreign     ',&
!!$         10x,'      self         foreign     ',      /,&
!!$             '    cm-1  ',&
!!$             '       1/(cm-1 molec/cm**2)    ',&
!!$         10x,'       1/(molec/cm**2) '        ,//)
!
      xkt=tave/radcn2

!!$      do 200 i=1,npth
!!$      vi=v1h+float(i-1)*dvh
!!$      if (vi.ge.v1abs .and. vi.le.v2abs) then
!!$         radfld=radfn(vi,xkt)
!!$         csh2or=csh2o(i) * radfld
!!$         cfh2or=cfh2o(i) * radfld
!!$         write (ipu,930) vi, csh2o(i), cfh2o(i), csh2or, cfh2or
!!$  930    format(f10.2, 1p, 2e15.4,10x, 1p, 2e15.4)
!!$      endif
!!$  200 continue

      do i = 1, npth
        vi=v1h+float(i-1)*dvh
        if (vi.ge.v1abs .and. vi.le.v2abs) then
          a_radfld(i) = radfn(vi,xkt)
          a_csh2o (i) = csh2o(i)
          a_cfh2o (i) = cfh2o(i)
        else
          a_radfld(i) = -1.0d0
          a_csh2o (i) = -1.0d0
          a_cfh2o (i) = -1.0d0
        end if
      end do

      if ( Arg_arraysize < npth ) then
        write( 6, * ) 'Array size is too small.', Arg_arraysize, npth
      end if
      Arg_NH2OArray = npth
      Arg_v1h       = v1h
      Arg_dvh       = dvh

    end subroutine call_mtckd

end module test06_mtckd_wrapper

Include 'cntnm_progr_procedures.f90'
