srtm_setcoef.F90

Go to the documentation of this file.

SUBROUTINE srtm_setcoef &
 & ( klev    , knmol ,&
 &   pavel   , ptavel   , pz      , ptz     , ptbound  ,&
 &   pcoldry , pwkl     ,&
 &   klaytrop, klayswtch, klaylow ,&
 &   pco2mult, pcolch4  , pcolco2 , pcolh2o , pcolmol  , pcoln2o  , pcolo2 , pcolo3 ,&
 &   pforfac , pforfrac , kindfor , pselffac, pselffrac, kindself ,&
 &   pfac00  , pfac01   , pfac10  , pfac11  ,&
 &   kjp     , kjt      , kjt1    &
 & )  

!     J. Delamere, AER, Inc. (version 2.5, 02/04/01)

!     Modifications:
!     JJMorcrette 030224   rewritten / adapted to ECMWF F90 system
!        M.Hamrud      01-Oct-2003 CY28 Cleaning

!     Purpose:  For a given atmosphere, calculate the indices and
!     fractions related to the pressure and temperature interpolations.

USE parkind1  ,ONLY : jpim     ,jprb
USE yomhook   ,ONLY : lhook,   dr_hook

USE parsrtm , ONLY : jplay
USE yoesrtwn, ONLY :  preflog, tref
!!  USE YOESWN  , ONLY : NDBUG
                   
IMPLICIT NONE

!-- Input arguments

INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV 
INTEGER(KIND=JPIM)               :: KNMOL ! Argument NOT used
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAVEL(jplay) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAVEL(jplay) 
REAL(KIND=JPRB)                  :: PZ(0:jplay) ! Argument NOT used
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTZ(0:jplay) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTBOUND 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCOLDRY(jplay) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PWKL(35,jplay) 
INTEGER(KIND=JPIM),INTENT(OUT)   :: KLAYTROP 
INTEGER(KIND=JPIM),INTENT(OUT)   :: KLAYSWTCH 
INTEGER(KIND=JPIM),INTENT(OUT)   :: KLAYLOW 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCO2MULT(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLCH4(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLCO2(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLH2O(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLMOL(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLN2O(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLO2(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLO3(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFORFAC(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFORFRAC(jplay) 
INTEGER(KIND=JPIM),INTENT(OUT)   :: KINDFOR(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSELFFAC(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSELFFRAC(jplay) 
INTEGER(KIND=JPIM),INTENT(OUT)   :: KINDSELF(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFAC00(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFAC01(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFAC10(jplay) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFAC11(jplay) 
INTEGER(KIND=JPIM),INTENT(OUT)   :: KJP(jplay) 
INTEGER(KIND=JPIM),INTENT(OUT)   :: KJT(jplay) 
INTEGER(KIND=JPIM),INTENT(OUT)   :: KJT1(jplay) 
!-- Output arguments

!-- local integers

INTEGER(KIND=JPIM) :: I_NLAYERS, INDBOUND, INDLEV0, JK
INTEGER(KIND=JPIM) :: JP1

!-- local reals

REAL(KIND=JPRB) :: Z_STPFAC, Z_TBNDFRAC, Z_T0FRAC, Z_PLOG, Z_FP, Z_FT, Z_FT1, Z_WATER, Z_SCALEFAC
REAL(KIND=JPRB) :: Z_FACTOR, Z_CO2REG, Z_COMPFP
REAL(KIND=JPRB) :: ZHOOK_HANDLE




IF (lhook) CALL dr_hook('SRTM_SETCOEF',0,zhook_handle)
i_nlayers = klev

z_stpfac = 296._jprb/1013._jprb

indbound = ptbound - 159._jprb
z_tbndfrac = ptbound - int(ptbound)
indlev0  = ptz(0) - 159._jprb
z_t0frac   = ptz(0) - int(ptz(0))

klaytrop  = 0
klayswtch = 0
klaylow   = 0

!IF (NDBUG.LE.3) THEN
!  print *,'-------- Computed in SETCOEF --------'
!  print 8990
8990 format(18x,'  T     PFAC00,    01,    10,    11  PCO2MULT     MOL   &
 & CH4      CO2      H2O      N2O      O2      O3      SFAC  &
 & SFRAC    FFAC    FFRAC  ISLF IFOR')  
!END IF

DO jk = 1, i_nlayers
!        Find the two reference pressures on either side of the
!        layer pressure.  Store them in JP and JP1.  Store in FP the
!        fraction of the difference (in ln(pressure)) between these
!        two values that the layer pressure lies.

  z_plog = log(pavel(jk))
  kjp(jk) = int(36. - 5*(z_plog+0.04))
  IF (kjp(jk) < 1) THEN
    kjp(jk) = 1
  ELSEIF (kjp(jk) > 58) THEN
    kjp(jk) = 58
  ENDIF
  jp1 = kjp(jk) + 1
  z_fp = 5. * (preflog(kjp(jk)) - z_plog)

!        Determine, for each reference pressure (JP and JP1), which
!        reference temperature (these are different for each  
!        reference pressure) is nearest the layer temperature but does
!        not exceed it.  Store these indices in JT and JT1, resp.
!        Store in FT (resp. FT1) the fraction of the way between JT
!        (JT1) and the next highest reference temperature that the 
!        layer temperature falls.

  kjt(jk) = int(3. + (ptavel(jk)-tref(kjp(jk)))/15.)
  IF (kjt(jk) < 1) THEN
    kjt(jk) = 1
  ELSEIF (kjt(jk) > 4) THEN
    kjt(jk) = 4
  ENDIF
  z_ft = ((ptavel(jk)-tref(kjp(jk)))/15.) - REAL(kjt(jk)-3)
  kjt1(jk) = int(3. + (ptavel(jk)-tref(jp1))/15.)
  IF (kjt1(jk) < 1) THEN
    kjt1(jk) = 1
  ELSEIF (kjt1(jk) > 4) THEN
    kjt1(jk) = 4
  ENDIF
  z_ft1 = ((ptavel(jk)-tref(jp1))/15.) - REAL(kjt1(jk)-3)

  z_water = pwkl(1,jk)/pcoldry(jk)
  z_scalefac = pavel(jk) * z_stpfac / ptavel(jk)

!        If the pressure is less than ~100mb, perform a different
!        set of species interpolations.

  IF (z_plog <= 4.56) GO TO 5300
  klaytrop =  klaytrop + 1
  IF (z_plog >= 6.62) klaylow = klaylow + 1

!        Set up factors needed to separately include the water vapor
!        foreign-continuum in the calculation of absorption coefficient.

  pforfac(jk) = z_scalefac / (1.+z_water)
  z_factor = (332.0-ptavel(jk))/36.0
  kindfor(jk) = min(2, max(1, int(z_factor)))
  pforfrac(jk) = z_factor - REAL(kindfor(jk))

!        Set up factors needed to separately include the water vapor
!        self-continuum in the calculation of absorption coefficient.

  pselffac(jk) = z_water * pforfac(jk)
  z_factor = (ptavel(jk)-188.0)/7.2
  kindself(jk) = min(9, max(1, int(z_factor)-7))
  pselffrac(jk) = z_factor - REAL(KINDSELF(JK) + 7)

!        Calculate needed column amounts.

  pcolh2o(jk) = 1.e-20 * pwkl(1,jk)
  pcolco2(jk) = 1.e-20 * pwkl(2,jk)
  pcolo3(jk) = 1.e-20 * pwkl(3,jk)
!         COLO3(LAY) = 0.
!         COLO3(LAY) = colo3(lay)/1.16
  pcoln2o(jk) = 1.e-20 * pwkl(4,jk)
  pcolch4(jk) = 1.e-20 * pwkl(6,jk)
  pcolo2(jk) = 1.e-20 * pwkl(7,jk)
  pcolmol(jk) = 1.e-20 * pcoldry(jk) + pcolh2o(jk)
!         colco2(lay) = 0.
!         colo3(lay) = 0.
!         coln2o(lay) = 0.
!         colch4(lay) = 0.
!         colo2(lay) = 0.
!         colmol(lay) = 0.
  IF (pcolco2(jk) == 0.) pcolco2(jk) = 1.e-32 * pcoldry(jk)
  IF (pcoln2o(jk) == 0.) pcoln2o(jk) = 1.e-32 * pcoldry(jk)
  IF (pcolch4(jk) == 0.) pcolch4(jk) = 1.e-32 * pcoldry(jk)
  IF (pcolo2(jk) == 0.) pcolo2(jk) = 1.e-32 * pcoldry(jk)
!        Using E = 1334.2 cm-1.
  z_co2reg = 3.55e-24 * pcoldry(jk)
  pco2mult(jk)= (pcolco2(jk) - z_co2reg) * &
   & 272.63*exp(-1919.4/ptavel(jk))/(8.7604e-4*ptavel(jk))  
  GO TO 5400

!        Above LAYTROP.
  5300    CONTINUE

!        Set up factors needed to separately include the water vapor
!        foreign-continuum in the calculation of absorption coefficient.

  pforfac(jk) = z_scalefac / (1.+z_water)
  z_factor = (ptavel(jk)-188.0)/36.0
  kindfor(jk) = 3
  pforfrac(jk) = z_factor - 1.0

!        Calculate needed column amounts.

  pcolh2o(jk) = 1.e-20 * pwkl(1,jk)
  pcolco2(jk) = 1.e-20 * pwkl(2,jk)
  pcolo3(jk)  = 1.e-20 * pwkl(3,jk)
  pcoln2o(jk) = 1.e-20 * pwkl(4,jk)
  pcolch4(jk) = 1.e-20 * pwkl(6,jk)
  pcolo2(jk)  = 1.e-20 * pwkl(7,jk)
  pcolmol(jk) = 1.e-20 * pcoldry(jk) + pcolh2o(jk)
  IF (pcolco2(jk) == 0.) pcolco2(jk) = 1.e-32 * pcoldry(jk)
  IF (pcoln2o(jk) == 0.) pcoln2o(jk) = 1.e-32 * pcoldry(jk)
  IF (pcolch4(jk) == 0.) pcolch4(jk) = 1.e-32 * pcoldry(jk)
  IF (pcolo2(jk) == 0.) pcolo2(jk)  = 1.e-32 * pcoldry(jk)
  z_co2reg = 3.55e-24 * pcoldry(jk)
  pco2mult(jk)= (pcolco2(jk) - z_co2reg) * &
   & 272.63*exp(-1919.4/ptavel(jk))/(8.7604e-4*ptavel(jk))  

  pselffac(jk) =0.0_jprb
  pselffrac(jk)=0.0_jprb
  kindself(jk) = 0

  5400    CONTINUE

!        We have now isolated the layer ln pressure and temperature,
!        between two reference pressures and two reference temperatures 
!        (for each reference pressure).  We multiply the pressure 
!        fraction FP with the appropriate temperature fractions to get 
!        the factors that will be needed for the interpolation that yields
!        the optical depths (performed in routines TAUGBn for band n).

  z_compfp = 1. - z_fp
  pfac10(jk) = z_compfp * z_ft
  pfac00(jk) = z_compfp * (1. - z_ft)
  pfac11(jk) = z_fp * z_ft1
  pfac01(jk) = z_fp * (1. - z_ft1)

!  IF (NDBUG.LE.3) THEN
!    print 9000,LAY,LAYTROP,JP(LAY),JT(LAY),JT1(LAY),TAVEL(LAY) &
!      &,FAC00(LAY),FAC01(LAY),FAC10(LAY),FAC11(LAY) &
!      &,CO2MULT(LAY),COLMOL(LAY),COLCH4(LAY),COLCO2(LAY),COLH2O(LAY) &
!      &,COLN2O(LAY),COLO2(LAY),COLO3(LAY),SELFFAC(LAY),SELFFRAC(LAY) &
!      &,FORFAC(LAY),FORFRAC(LAY),INDSELF(LAY),INDFOR(LAY)
  9000 format(1x,2i3,3i4,f6.1,4f7.2,12e9.2,2i5)
!  END IF

ENDDO

!----------------------------------------------------------------------- 
IF (lhook) CALL dr_hook('SRTM_SETCOEF',1,zhook_handle)
END SUBROUTINE srtm_setcoef