radae.f90

来自「CCSM Research Tools: Community Atmospher」· F90 代码 · 共 1,421 行 · 第 1/5 页

               emm(i,2)    = 0.5*(co2em(i,k2) + co2eml(i,k2))
               tbar(i,3)   = 0.5*(tbar(i,2) + tbar(i,1))
               emm(i,3)    = emm(i,1)
               tbar(i,4)   = tbar(i,3)
               emm(i,4)    = emm(i,2)
               o3emm(i,1)  = 0.5*(dbvtit(i,k2+1) + dbvtly(i,k2))
               o3emm(i,2)  = 0.5*(dbvtit(i,k2) + dbvtly(i,k2))
               o3emm(i,3)  = o3emm(i,1)
               o3emm(i,4)  = o3emm(i,2)
               temh2o(i,1) = tbar(i,1)
               temh2o(i,2) = tbar(i,2)
               temh2o(i,3) = tbar(i,1)
               temh2o(i,4) = tbar(i,2)
               dpnm(i)     = pnm(i,k2+1) - pnm(i,k2)
            end do
!
!  Weighted Planck functions for trace gases
!
            do wvl = 1,14
               do i = 1,ncol
                  bplnk(wvl,i,1) = 0.5*(abplnk1(wvl,i,k2+1) + abplnk2(wvl,i,k2))
                  bplnk(wvl,i,2) = 0.5*(abplnk1(wvl,i,k2) + abplnk2(wvl,i,k2))
                  bplnk(wvl,i,3) = bplnk(wvl,i,1)
                  bplnk(wvl,i,4) = bplnk(wvl,i,2)
               end do
            end do

            do i=1,ncol
               rdpnmsq    = 1./(pnmsq(i,k2+1) - pnmsq(i,k2))
               rdpnm      = 1./dpnm(i)
               p1         = .5*(pbr(i,k2) + pnm(i,k2+1))
               p2         = .5*(pbr(i,k2) + pnm(i,k2  ))
               uinpl(i,1) =  (pnmsq(i,k2+1) - p1**2)*rdpnmsq
               uinpl(i,2) = -(pnmsq(i,k2  ) - p2**2)*rdpnmsq
               uinpl(i,3) = -(pnmsq(i,k2  ) - p1**2)*rdpnmsq
               uinpl(i,4) =  (pnmsq(i,k2+1) - p2**2)*rdpnmsq
               winpl(i,1) = (.5*( pnm(i,k2+1) - pbr(i,k2)))*rdpnm
               winpl(i,2) = (.5*(-pnm(i,k2  ) + pbr(i,k2)))*rdpnm
               winpl(i,3) = (.5*( pnm(i,k2+1) + pbr(i,k2)) - pnm(i,k2  ))*rdpnm
               winpl(i,4) = (.5*(-pnm(i,k2  ) - pbr(i,k2)) + pnm(i,k2+1))*rdpnm
               tmp1       = 1./(piln(i,k2+1) - piln(i,k2))
               tmp2       = piln(i,k2+1) - pmln(i,k2)
               tmp3       = piln(i,k2  ) - pmln(i,k2)
               zinpl(i,1) = (.5*tmp2          )*tmp1
               zinpl(i,2) = (        - .5*tmp3)*tmp1
               zinpl(i,3) = (.5*tmp2 -    tmp3)*tmp1
               zinpl(i,4) = (   tmp2 - .5*tmp3)*tmp1
               pinpl(i,1) = 0.5*(p1 + pnm(i,k2+1))
               pinpl(i,2) = 0.5*(p2 + pnm(i,k2  ))
               pinpl(i,3) = 0.5*(p1 + pnm(i,k2  ))
               pinpl(i,4) = 0.5*(p2 + pnm(i,k2+1))
            end do
            do 400 kn=1,4
               do i=1,ncol
                  u(i)     = uinpl(i,kn)*abs(plh2o(i,k2) - plh2o(i,k2+1))
                  sqrtu(i) = sqrt(u(i))
                  dw(i)    = abs(w(i,k2) - w(i,k2+1))
                  pnew(i)  = u(i)/(winpl(i,kn)*dw(i))
                  pnew_mks  = pnew(i) * sslp_mks
                  t_p = min(max(tbar(i,kn), min_tp_h2o), max_tp_h2o)
                  iest = floor(t_p) - min_tp_h2o
                  esx = estblh2o(iest) + (estblh2o(iest+1)-estblh2o(iest)) * &
                       (t_p - min_tp_h2o - iest)
                  qsx = epsilo * esx / (pnew_mks - omeps * esx)
                  q_path = dw(i) / ABS(dpnm(i)) / rga

                  ds2c     = abs(s2c(i,k2) - s2c(i,k2+1))
                  uc1(i)   = uinpl(i,kn)*ds2c
                  pch2o    = uc1(i)
                  uc1(i)   = (uc1(i) + 1.7e-3*u(i))*(1. +  2.*uc1(i))/(1. + 15.*uc1(i))
                  dtx(i)      = temh2o(i,kn) - 250.
                  dty(i)      = tbar(i,kn) - 250.

                  fwk(i)    = fwcoef + fwc1/(1. + fwc2*u(i))
                  fwku(i)   = fwk(i)*u(i)
!
! Define variables for C/H/E fit
!
! abso(i,1)     0 -  800 cm-1   h2o rotation band
! abso(i,1)  1200 - 2200 cm-1   h2o vibration-rotation band
! abso(i,2)   800 - 1200 cm-1   h2o window
!
! Separation between rotation and vibration-rotation dropped, so
!                only 2 slots needed for H2O absorptivity
!
! Notation:
! U   = integral (P/P_0 dW)  
! P   = atmospheric pressure
! P_0 = reference atmospheric pressure
! W   = precipitable water path
! T_e = emission temperature
! T_p = path temperature
! RH  = path relative humidity
!
!
! Terms for asymptotic value of emissivity
!
                  te1  = temh2o(i,kn)
                  te2  = te1 * te1
                  te3  = te2 * te1
                  te4  = te3 * te1
                  te5  = te4 * te1

!
! Indices for lines and continuum tables 
! Note: because we are dealing with the nearest layer,
!       the Hulst-Curtis-Godson corrections
!       for inhomogeneous paths are not applied.
!
                  uvar = u(i)*fdif
                  log_u  = min(log10(max(uvar, min_u_h2o)), max_lu_h2o)
                  dvar = (log_u - min_lu_h2o) / dlu_h2o
                  iu = min(max(int(aint(dvar,r8)) + 1, 1), n_u - 1)
                  iu1 = iu + 1
                  wu = dvar - floor(dvar)
                  wu1 = 1.0 - wu
                  
                  log_p  = min(log10(max(pnew(i), min_p_h2o)), max_lp_h2o)
                  dvar = (log_p - min_lp_h2o) / dlp_h2o
                  ip = min(max(int(aint(dvar,r8)) + 1, 1), n_p - 1)
                  ip1 = ip + 1
                  wp = dvar - floor(dvar)
                  wp1 = 1.0 - wp
                  
                  dvar = (t_p - min_tp_h2o) / dtp_h2o
                  itp = min(max(int(aint(dvar,r8)) + 1, 1), n_tp - 1)
                  itp1 = itp + 1
                  wtp = dvar - floor(dvar)
                  wtp1 = 1.0 - wtp
                  
                  t_e = min(max(temh2o(i,kn)-t_p,min_te_h2o),max_te_h2o)
                  dvar = (t_e - min_te_h2o) / dte_h2o
                  ite = min(max(int(aint(dvar,r8)) + 1, 1), n_te - 1)
                  ite1 = ite + 1
                  wte = dvar - floor(dvar)
                  wte1 = 1.0 - wte
                  
                  rh_path = min(max(q_path / qsx, min_rh_h2o), max_rh_h2o)
                  dvar = (rh_path - min_rh_h2o) / drh_h2o
                  irh = min(max(int(aint(dvar,r8)) + 1, 1), n_rh - 1)
                  irh1 = irh + 1
                  wrh = dvar - floor(dvar)
                  wrh1 = 1.0 - wrh

                  w_0_0_ = wtp  * wte
                  w_0_1_ = wtp  * wte1
                  w_1_0_ = wtp1 * wte 
                  w_1_1_ = wtp1 * wte1
                  
                  w_0_00 = w_0_0_ * wrh
                  w_0_01 = w_0_0_ * wrh1
                  w_0_10 = w_0_1_ * wrh
                  w_0_11 = w_0_1_ * wrh1
                  w_1_00 = w_1_0_ * wrh
                  w_1_01 = w_1_0_ * wrh1
                  w_1_10 = w_1_1_ * wrh
                  w_1_11 = w_1_1_ * wrh1
                        
                  w00_00 = wp  * w_0_00 
                  w00_01 = wp  * w_0_01 
                  w00_10 = wp  * w_0_10 
                  w00_11 = wp  * w_0_11 
                  w01_00 = wp  * w_1_00 
                  w01_01 = wp  * w_1_01 
                  w01_10 = wp  * w_1_10 
                  w01_11 = wp  * w_1_11 
                  w10_00 = wp1 * w_0_00 
                  w10_01 = wp1 * w_0_01 
                  w10_10 = wp1 * w_0_10 
                  w10_11 = wp1 * w_0_11 
                  w11_00 = wp1 * w_1_00 
                  w11_01 = wp1 * w_1_01 
                  w11_10 = wp1 * w_1_10 
                  w11_11 = wp1 * w_1_11 

!
! Non-window absorptivity
!
                  ib = 1

                  fa = fat(1,ib) + &
                       fat(2,ib) * te1 + &
                       fat(3,ib) * te2 + &
                       fat(4,ib) * te3 + &
                       fat(5,ib) * te4 + &
                       fat(6,ib) * te5
      
                  a_star = &
                       ah2onw(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &
                       ah2onw(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &
                       ah2onw(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &
                       ah2onw(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &
                       ah2onw(ip , itp , iu1, ite , irh ) * w11_11 * wu  + &
                       ah2onw(ip , itp , iu1, ite , irh1) * w11_10 * wu  + &
                       ah2onw(ip , itp , iu1, ite1, irh ) * w11_01 * wu  + &
                       ah2onw(ip , itp , iu1, ite1, irh1) * w11_00 * wu  + &
                       ah2onw(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &
                       ah2onw(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &
                       ah2onw(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &
                       ah2onw(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &
                       ah2onw(ip , itp1, iu1, ite , irh ) * w10_11 * wu  + &
                       ah2onw(ip , itp1, iu1, ite , irh1) * w10_10 * wu  + &
                       ah2onw(ip , itp1, iu1, ite1, irh ) * w10_01 * wu  + &
                       ah2onw(ip , itp1, iu1, ite1, irh1) * w10_00 * wu  + &
                       ah2onw(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &
                       ah2onw(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &
                       ah2onw(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &
                       ah2onw(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &
                       ah2onw(ip1, itp , iu1, ite , irh ) * w01_11 * wu  + &
                       ah2onw(ip1, itp , iu1, ite , irh1) * w01_10 * wu  + &
                       ah2onw(ip1, itp , iu1, ite1, irh ) * w01_01 * wu  + &
                       ah2onw(ip1, itp , iu1, ite1, irh1) * w01_00 * wu  + &
                       ah2onw(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &
                       ah2onw(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &
                       ah2onw(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &
                       ah2onw(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &
                       ah2onw(ip1, itp1, iu1, ite , irh ) * w00_11 * wu  + &
                       ah2onw(ip1, itp1, iu1, ite , irh1) * w00_10 * wu  + &
                       ah2onw(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu  + &
                       ah2onw(ip1, itp1, iu1, ite1, irh1) * w00_00 * wu
 
                  abso(i,ib) = min(max(fa * a_star, 0.0_r8), 1.0_r8)
!
! Invoke linear limit for scaling wrt u below min_u_h2o
!
                  if (uvar < min_u_h2o) then
                     uscl = uvar / min_u_h2o
                     abso(i,ib) = abso(i,ib) * uscl
                  endif
                               
!
! Window absorptivity
!
                  ib = 2

                  fa = fat(1,ib) + &
                       fat(2,ib) * te1 + &
                       fat(3,ib) * te2 + &
                       fat(4,ib) * te3 + &
                       fat(5,ib) * te4 + &
                       fat(6,ib) * te5
      
                  a_star = &
                       ah2ow(ip , itp , iu , ite , irh ) * w11_11 * wu1 + &
                       ah2ow(ip , itp , iu , ite , irh1) * w11_10 * wu1 + &
                       ah2ow(ip , itp , iu , ite1, irh ) * w11_01 * wu1 + &
                       ah2ow(ip , itp , iu , ite1, irh1) * w11_00 * wu1 + &
                       ah2ow(ip , itp , iu1, ite , irh ) * w11_11 * wu  + &
                       ah2ow(ip , itp , iu1, ite , irh1) * w11_10 * wu  + &
                       ah2ow(ip , itp , iu1, ite1, irh ) * w11_01 * wu  + &
                       ah2ow(ip , itp , iu1, ite1, irh1) * w11_00 * wu  + &
                       ah2ow(ip , itp1, iu , ite , irh ) * w10_11 * wu1 + &
                       ah2ow(ip , itp1, iu , ite , irh1) * w10_10 * wu1 + &
                       ah2ow(ip , itp1, iu , ite1, irh ) * w10_01 * wu1 + &
                       ah2ow(ip , itp1, iu , ite1, irh1) * w10_00 * wu1 + &
                       ah2ow(ip , itp1, iu1, ite , irh ) * w10_11 * wu  + &
                       ah2ow(ip , itp1, iu1, ite , irh1) * w10_10 * wu  + &
                       ah2ow(ip , itp1, iu1, ite1, irh ) * w10_01 * wu  + &
                       ah2ow(ip , itp1, iu1, ite1, irh1) * w10_00 * wu  + &
                       ah2ow(ip1, itp , iu , ite , irh ) * w01_11 * wu1 + &
                       ah2ow(ip1, itp , iu , ite , irh1) * w01_10 * wu1 + &
                       ah2ow(ip1, itp , iu , ite1, irh ) * w01_01 * wu1 + &
                       ah2ow(ip1, itp , iu , ite1, irh1) * w01_00 * wu1 + &
                       ah2ow(ip1, itp , iu1, ite , irh ) * w01_11 * wu  + &
                       ah2ow(ip1, itp , iu1, ite , irh1) * w01_10 * wu  + &
                       ah2ow(ip1, itp , iu1, ite1, irh ) * w01_01 * wu  + &
                       ah2ow(ip1, itp , iu1, ite1, irh1) * w01_00 * wu  + &
                       ah2ow(ip1, itp1, iu , ite , irh ) * w00_11 * wu1 + &
                       ah2ow(ip1, itp1, iu , ite , irh1) * w00_10 * wu1 + &
                       ah2ow(ip1, itp1, iu , ite1, irh ) * w00_01 * wu1 + &
                       ah2ow(ip1, itp1, iu , ite1, irh1) * w00_00 * wu1 + &
                       ah2ow(ip1, itp1, iu1, ite , irh ) * w00_11 * wu  + &
                       ah2ow(ip1, itp1, iu1, ite , irh1) * w00_10 * wu  + &
                       ah2ow(ip1, itp1, iu1, ite1, irh ) * w00_01 * wu  + &
                       ah2ow(ip1, itp1, iu1, ite1, irh1) * w00_00 * wu
 
                  abso(i,ib) = min(max(fa * a_star, 0.0_r8), 1.0_r8) 
!
! Invoke linear limit for scaling wrt u below min_u_h2o
!
                  if (uvar < min_u_