!= MT_CKD continuum model サブルーチン化プログラム
!
!= subroutine from MT_CKD continuum model 
!
! サブルーチン化: 大西将徳
! 2014/05/29
! 2014/06/04
! 2016/07/21 bug fix 

!     path:      $Source: /project/rc/rc1/cvsroot/rc/cntnm/src/cntnm_progr.f,v $
!     author:    $Author: jdelamer $
!     revision:  $Revision: 1.2 $
!     created:   $Date: 2011/03/29 20:16:57 $
!
!
!  --------------------------------------------------------------------------
! |  Copyright ©, Atmospheric and Environmental Research, Inc., 2011         |
! |                                                                          |
! |  All rights reserved. This source code is part of the MT_CKD continuum   |
! |  software and is designed for scientific and research purposes.          |
! |  Atmospheric and Environmental Research, Inc. (AER) grants USER          |
! |  the right to download, install, use and copy this software              |
! |  for scientific and research purposes only. This software may be         |
! |  redistributed as long as this copyright notice is reproduced on any     |
! |  copy made and appropriate acknowledgment is given to AER. This          |
! |  software or any modified version of this software may not be            |
! |  incorporated into proprietary software or commercial software           |
! |  offered for sale.                                                       |
! |                                                                          |
! |  This software is provided as is without any express or implied          |
! |  warranties.                                                             |
! |                       (http://www.rtweb.aer.com/)                        |
!  --------------------------------------------------------------------------
!
!                  PROGRAM DRCNTNM
  subroutine MT_CKD (InPress, InTemp, ArrVMR, WNmin, WNmax, WNres, AbsorpCS)
!
!
!
!     The mt_ckd water vapor continuum is a completely new continuum  
!     formulation based on a collision induced  component and a sub-Lorentzian
!     line wing component.  Both the water vapor continuum coefficients and
!     those for other molecules are constrained to agree with accurate
!     measurements of continuum absorption in the spectral regions where such
!     measurements exist.
!
!     This is an updated version of the continuum program:
!     this version provides optical depths on file CNTNM.OPTDPT as before:
!     it also provides the continuum coefficients on file  WATER.COEF
!
!     the length of the header records may vary by version:
!         in this version the WATER.COEF header information is 47 records
!         in this version the CNTNM.OPTDT header information is 34 records
!
!     presumably the user will want to create an input file to address
!     individual requirements
!
!
      IMPLICIT REAL*8           (V)
!
      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)
!
      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/

      !================= 2014/06/02 onishi
      real(8) :: InPress
      real(8) :: InTemp
      real(8) :: ArrVMR(1:5) !H2O, CO2, O3, O2, N2
      real(8) :: WNmin
      real(8) :: WNmax
      real(8) :: WNres
      real(8) :: AbsorpCS(1:2, 1:5001)
      !=================- 2013/03/07 onishi
      real(8) :: vmrh2o  !2014/06/02 onishi
      real(8) :: vmrco2
      real(8) :: vmro3
      real(8) :: vmrn2   !2014/01/15 onishi
      real(8) :: vmrar   !2014/01/15 onishi
      real(8) :: vmro2   !2014/01/15 onishi
      !real(8) :: USvalue(0:9) ! 2013/07/12 onishi 2014/01/15 onishi rev
      !integer :: ifo = 11  !! SET
      !integer :: fot = 12  !2014/05/29
      ! ilay =        45 !!
      !=================  2013/03/07 onishi
 !     open(fot, file='test_out')
!      open(ifo, FILE='InputMTCKD')
!      do i=1, 1
!        read(ifo, *)
!      end do

!      do ilay = 1, 10 ! 2013/11/18 onishi
!      read(ifo, *) USvalue
!      print *, USvalue
!      !================== 2013/03/07 onishi
!
      RADCN1=2.*PLANCK*CLIGHT*CLIGHT*1.E-07
      RADCN2=PLANCK*CLIGHT/BOLTZ
!
      icflg = -999
!
      do 1, i=1,7
         xcnt(i)=1.
 1    continue

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

      do 3, i=1,60
         wk(i)=0.
 3    continue
!
      ird = 55
      ipr = 66
      ipu = 7
!                                SET
!      OPEN (ipr,FILE='CNTNM_US_layAll.OPTDPT')
!      OPEN (ipu,FILE='WATER.COEF')
!
      !print *
      !print *, '  This version is limited to 5000 values  '
!
!   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 = 1013.
      TAVE =  296.
!
      VMRH2O = 0.01
!
      xlength = 1.
!
      !print *
      !print *,' *** For this program, vmr_h2o is taken ', &
      ! 'with respect to the total column ***'

      !print *
      !print *,' read: pressure (mb)  if negative use default values'
!      read *, press_rd
!      press_rd = USvalue(1)
      press_rd = InPress

!!$      if (press_rd .gt. 0.) then
      if (press_rd .ge. 0.) then
         pave = press_rd
!         print *,' read:   temperature (K)'
!         read *, tave
         !tave = USvalue(2)
         tave = InTemp
!         print *,' read:   path length (cm)'
!         read *, xlength
         xlength = 100.0 !USvalue(3)
!         xlength = 10000.0 !USvalue(3)
!         print *,' read:   vmr h2o '
!         read *, vmrh2o
         vmrh2o = ArrVMR(1) !USvalue(4)
!         print *, ' read:  vmr CO2 '
!         read *, vmrco2
         vmrco2 = ArrVMR(2) !USvalue(5)
!         print *, ' read: vmr O3 '
!         read *, vmro3
         vmro3  = ArrVMR(3) !USvalue(6)

         vmro2  = ArrVMR(4) !USvalue(7) ! 2014/01/15 onishi
         vmrn2  = ArrVMR(5) !USvalue(8) ! 2014/01/15 onishi
         vmrar  = 0.0 ! 2014/06/02 USvalue(9) ! 2014/01/15 onishi
      endif
      !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)

!! == SET
!  pave = 955.9    !! pressure
!  tave = 284.95   !! tempereture
!  xlength = 1.0E+5!! vertical length [cm]
!  vmrh2o =7.69E-3 !! volume mixing ratio of H2O
!! == SET
!
!     It may be preferable to specifiy the water column directly!
!
      WTOT = XLOSMT*(PAVE/1013.)*(273./TAVE)* xlength
!      print *, WTOT ! onishi 2014/06/02
!      print *, xlength
!
!      W_dry = WTOT * (1.-VMRH2O)  
      W_dry = WTOT / (1.0 + VMRH2O)  ! onishi 2014/06/02 
!
!     ww is column of dry air;  vol. mix. ratios are based on dry air
!
! argon:
!      WA     = 0.009     * W_dry
      WA     = vmrar     * W_dry  ! 2014/01/15 onishi
! nitrogen:
!      WN2    = 0.78      * W_dry
      WN2    = vmrn2     * W_dry  ! 2014/01/15 onishi
! oxygen:
!      WK(7)  = 0.21      * W_dry
      WK(7)  = vmro2     * W_dry  ! 2014/01/15 onishi
! carbon dioxide:
      WK(2)  = vmrco2  * W_dry

!      WK(2) = 0.

! ozone ! 2013/01/12 onishi adds !! SET
      WK(3)  = 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
!                                !! SET ! 2013/03/07 onishi modify
      V1ABS = WNmin !0.0_8              !! 0.
      V2ABS = WNmax !20000.0_8 !26000.0_8          !! 10000.

      DVABS = WNres !20000.0_8/5000.0_8 !26000.0_8/5000.0_8    !! 2.
! ..........................................................
!      write (*,*) '  v1abs,  v2abs,  dvabs  '
!      read  (*,*)    v1abs,  v2abs,  dvabs
! ..........................................................

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

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

!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 (fot,980) (WK(M),M=1,7),WBROAD !20140529 onishi
!
      !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,//)
!
      jrad=1
!
      v1 = v1abs
      v2 = v2abs
!
      CALL CONTNM(JRAD)
!
!      DO 100 I=1,NPTABS
!      VI=V1ABS+FLOAT(I-1)*DVABS
!100   WRITE (ipr, 910) VI, ABSRB(I)

      do I=1,NPTABS
        VI=V1ABS+FLOAT(I-1)*DVABS
        !WRITE (ipr, 910) VI, ABSRB(I)
!        write (fot, 910) VI, ABSRB(I)   ! 20140529 onishi
        AbsorpCS(1,I) = VI              ! 20140603 onishi
        AbsorpCS(2,I) = (ABSRB(I)*1d0/WTOT) ![cm2/molecule]! 20140603 onishi !20160721 (*1d0)
        !print *, AbsorpCS(1:2,I)
        !print *, ABSRB(I)*1d0, WTOT*1d0, ABSRB(I)*1d0/WTOT
      end do
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

!      end do ! ilay ! 2013/07/18 onishi
!
      END SUBROUTINE
!**********************************************************************
      Block Data phys_consts
!
      COMMON /CONSTS/ PI,PLANCK,BOLTZ,CLIGHT,AVOGAD,ALOSMT,GASCON, &
                      RADCN1,RADCN2,GRAV,CPDAIR,AIRMWT,SECDY
!
!!      DATA PI /3.1415926535898 /
!
!    Constants from NIST 01/11/2002
!
!!      DATA PLANCK / 6.62606876E-27 /, BOLTZ  / 1.3806503E-16 /, &
!!           CLIGHT / 2.99792458E+10 /, &
!!           AVOGAD / 6.02214199E+23 /, ALOSMT / 2.6867775E+19 /, &
!!           GASCON / 8.314472E+07 / &
!!           RADCN1 / 1.191042722E-12 /, RADCN2 / 1.4387752    /


      DATA PI /3.1415926535898 /, &
           PLANCK / 6.62606876E-27 /, BOLTZ  / 1.3806503E-16 /, &
           CLIGHT / 2.99792458E+10 /, &
           AVOGAD / 6.02214199E+23 /, ALOSMT / 2.6867775E+19 /, &
           GASCON / 8.314472E+07 / &
           RADCN1 / 1.191042722E-12 /, RADCN2 / 1.4387752    /, &
           GRAV /9.8 /, CPDAIR / 0.0 /, & !20140529 onishi
           AIRMWT / 0.0 /, SECDY / 0.0/ ! 20140529 onishi
!
!     Pi was obtained from   PI = 2.*ASIN(1.)                             A03980
!
!     units are generally cgs
!
!     The first and second radiation constants are taken from NIST.
!     They were previously obtained from the relations:
!                            RADCN1 = 2.*PLANCK*CLIGHT*CLIGHT*1.E-07      A03990
!                            RADCN2 = PLANCK*CLIGHT/BOLTZ                 A04000
      end
!
      BLOCK DATA
!
      COMMON /LAMCHN/ ONEPL,ONEMI,EXPMIN,ARGMIN
!
      DATA ONEPL/1.001/, ONEMI/0.999/, ARGMIN/34./
!                                                                         A07710

      END
      BLOCK DATA cntnm
!
      IMPLICIT REAL*8           (V)
!
      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
!                                                                         F00100
      CHARACTER*8      XID,       HMOLID,      YID
      REAL*8               SECANT,       XALTZ
!
      character*8 holn2
!
      DATA HMOlid/ '  H2O   ' , '  CO2   ' , '   O3   ' , '  N2O   ' , &
                   '   CO   ' , '  CH4   ' , '   O2   ' , 53*'        '/
!
      DATA HOLN2 / ' OTHER'/
!
      end
      SUBROUTINE XINT (V1A,V2A,DVA,A,AFACT,VFT,DVR3,R3,N1R3,N2R3)
!                                                                         B17530
!                                                                         B17870
      IMPLICIT REAL*8           (V)
!                                                                         B17550
!     THIS SUBROUTINE INTERPOLATES THE A ARRAY STORED                     B17560
!     FROM V1A TO V2A IN INCREMENTS OF DVA USING A MULTIPLICATIVE         B17570
!     FACTOR AFACT, INTO THE R3 ARRAY FROM LOCATION N1R3 TO N2R3 IN       B17580
!     INCREMENTS OF DVR3                                                  B17590
!                                                                         B17600
      COMMON /LAMCHN/ ONEPL,ONEMI,EXPMIN,ARGMIN
      DIMENSION A(*),R3(*)
!                                                                         B17630
      RECDVA = 1./DVA
      ILO = (V1A+DVA-VFT)/DVR3+1.+ONEMI
      ILO = MAX(ILO,N1R3)
      IHI = (V2A-DVA-VFT)/DVR3+ONEMI
      IHI = MIN(IHI,N2R3)
!                                                                         B17690
      DO 10 I = ILO, IHI
         VI = VFT+DVR3*FLOAT(I-1)
         J = (VI-V1A)*RECDVA+ONEPL
         VJ = V1A+DVA*FLOAT(J-1)
         P = RECDVA*(VI-VJ)
         C = (3.-2.*P)*P*P
         B = 0.5*P*(1.-P)
         B1 = B*(1.-P)
         B2 = B*P
         CONTI = -A(J-1)*B1+A(J)*(1.-C+B2)+A(J+1)*(C+B1)-A(J+2)*B2
         R3(I) = R3(I)+CONTI*AFACT
   10 CONTINUE
!                                                                         B17820
      RETURN
!                                                                         B17840
      END
      FUNCTION RADFN (VI,XKT)
!                                                                         B17870
      IMPLICIT REAL*8           (V)
!                                                                         B17890
!     FUNCTION RADFN CALCULATES THE RADIATION TERM FOR THE LINE SHAPE     B17900
!                                                                         B17910
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   B17920
!                                                                         B17930
!               LAST MODIFICATION:    12 AUGUST 1991                      B17940
!                                                                         B17950
!                  IMPLEMENTATION:    R.D. WORSHAM                        B17960
!                                                                         B17970
!             ALGORITHM REVISIONS:    S.A. CLOUGH                         B17980
!                                     R.D. WORSHAM                        B17990
!                                     J.L. MONCET                         B18000
!                                                                         B18010
!                                                                         B18020
!                     ATMOSPHERIC AND ENVIRONMENTAL RESEARCH INC.         B18030
!                     840 MEMORIAL DRIVE,  CAMBRIDGE, MA   02139          B18040
!                                                                         B18050
!----------------------------------------------------------------------   B18060
!                                                                         B18070
!               WORK SUPPORTED BY:    THE ARM PROGRAM                     B18080
!                                     OFFICE OF ENERGY RESEARCH           B18090
!                                     DEPARTMENT OF ENERGY                B18100
!                                                                         B18110
!                                                                         B18120
!      SOURCE OF ORIGINAL ROUTINE:    AFGL LINE-BY-LINE MODEL             B18130
!                                                                         B18140
!                                             FASCOD3                     B18150
!                                                                         B18160
!                                                                         B18170
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   B18180
!                                                                         B18190
      COMMON /LAMCHN/ ONEPL,ONEMI,EXPMIN,ARGMIN
!                                                                         B18210
!      IN THE SMALL XVIOKT REGION 0.5 IS REQUIRED                         B18220
!                                                                         B18230
      XVI = VI
!                                                                         B18250
      IF (XKT.GT.0.0) THEN
!                                                                         B18270
         XVIOKT = XVI/XKT
!                                                                         B18290
         IF (XVIOKT.LE.0.01) THEN
            RADFN = 0.5*XVIOKT*XVI
!                                                                         B18320
         ELSEIF (XVIOKT.LE.10.0) THEN
            EXPVKT = EXP(-XVIOKT)
            RADFN = XVI*(1.-EXPVKT)/(1.+EXPVKT)
!                                                                         B18360
         ELSE
            RADFN = XVI
         ENDIF
!                                                                         B18400
      ELSE
         RADFN = XVI
      ENDIF
!                                                                         B18440
      RETURN
!                                                                         B18460
      END
!*******
!*******
!*******
      Include 'contnm.f90'
