!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module mcm_lw_mod ! TK modified from ajb code ! /net/ajb/radiation_code/updates/v4.3/lw_mod.F ! Added interface routine (lw_rad_ss) which is called by ! fsrad and which calls lwcool in this ! module after constructing the appropriate inputs. USE Constants_Mod, ONLY: grav, tfreeze Use Fms_Mod, ONLY: Error_Mesg, FATAL, & write_version_number, mpp_pe, mpp_root_pe implicit none private integer, parameter :: nb_lw=19, ng=3, ngp=ng+1 integer :: ix, jx, kx, kp, km !------------ VERSION NUMBER ---------------- character(len=128) :: version = '$Id: mcm_lw.F90,v 13.0 2006/03/28 21:09:36 fms Exp $' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' logical :: module_is_initialized = .false. public :: MCM_LW_RAD, mcm_lw_init, mcm_lw_end ! ------------------------------------------------- ! TK NOTE: not ready for this yet... implicit none !----------------------------------------------------------------------- !--------------------- G L O B A L D A T A --------------------------- !----------------------------------------------------------------------- ! #include "parm.h" ! TK note: (ix,jx are set up below -- depend on domain decomposition) ! integer, parameter :: ix = 96 ! integer, parameter :: jx = 80 ! TK note: (kx,kp,km defined in rdparm) ! integer, parameter :: kx = 14 ! integer, parameter :: kp = kx + 1 ! integer, parameter :: km = kx - 1 ! TK note: (nb,ng,ngp are defined in rdparm) ! integer, parameter :: nb=19 ! integer, parameter :: ng=3 ! integer, parameter :: ngp=ng +1 ! nb=number of spectral bands ! ng=number of absorbing gases ! TK: This is a name change to avoid having to change nb locally: INTEGER, PARAMETER :: nb = nb_lw INTEGER :: LMAX, i real :: grav_accel real :: pi_alpha real :: co2_mixrat real :: t_freeze real :: h2o_line_width(nb) real :: h2o_line_strength(nb) real :: h2o_corr_a(nb,2) real :: h2o_corr_b(nb,2) real :: absorp_table(30,12) real :: d_absorp_dt(7,12) data h2o_line_width/ & .93000e-01, .18200e-00, .94000e-01, .79700e-01, .73300e-01, & .52000e-01, .67000e-01, .45900e-01, .10000e+01, .10000e+01, & .89000e-01, .23000e-00, .32000e-00, .29600e-00, .45200e-00, & .35900e+00, .16500e+00, .10400e+00, .11600e+00 / data h2o_line_strength/ & .57975e+03, .72103e+04, .60248e+04, .14039e+04, .78952e+02, & .42700e+01, .71500e-01, .27500e-01, .00000e+00, .00000e+00, & .12650e+02, .13440e+03, .63290e+03, .33120e+03, .43410e+03, & .13600e+03, .35650e+02, .90150e+01, .15290e+01 / data (h2o_corr_a(i,1),i=1,19)/ & -.87600e-02, -.26800e-02, .20300e-02, .96100e-02, .15820e-01, & .01705e+00, .02620e+00, .02918e+00, .0 , .0 , & .0 , .0 , .0 , .0 , .0 , & .0 , .0 , .0 , .0 / data (h2o_corr_a(i,2),i=1,19)/ & -.67500e-02, -.29300e-02, .14300e-02, .98400e-02, .13710e-01, & .01579e+00, .02410e+00, .02596e+00, .0 , .0 , & .0 , .0 , .0 , .0 , .0 , & .0 , .0 , .0 , .0 / data (h2o_corr_b(i,1),i=1,19)/ & .14150e-04, .15700e-05, -.10300e-04, -.43140e-04, -.37440e-04, & -.51440e-04, -.74100e-04, -.80760e-04, .0 , .0 , & .0 , .0 , .0 , .0 , .0 , & .0 , .0 , .0 , .0 / data (h2o_corr_b(i,2),i=1,19)/ & .85500e-05, .20100e-05, -.13000e-04, -.40810e-04, -.16150e-04, & -.44510e-04, -.40300e-04, -.66720e-04, .0 , .0 , & .0 , .0 , .0 , .0 , .0 , & .0 , .0 , .0 , .0 / data (absorp_table(i,12),i=1,30)/ & .0000,.0001,.0001,.0002,.0004,.0006,.0011,.0020,.0034,.0058, & .0097,.0156,.0248,.0386,.0586,.0861,.1205,.1607,.2040,.2480, & .2918,.3362,.3812,.4250,.4665,.5045,.5386,.5707,.6037,.6395/ data (absorp_table(i,11),i=1,30)/ & .0000,.0001,.0001,.0002,.0004,.0006,.0011,.0020,.0034,.0057, & .0093,.0149,.0231,.0351,.0520,.0744,.1032,.1388,.1803,.2255, & .2715,.3157,.3580,.3993,.4404,.4812,.5200,.5558,.5899,.6248/ data (absorp_table(i,10),i=1,30)/ & .0000,.0001,.0001,.0002,.0004,.0006,.0011,.0019,.0033,.0054, & .0086,.0131,.0194,.0280,.0398,.0557,.0771,.1047,.1385,.1777, & .2206,.2648,.3080,.3497,.3910,.4328,.4748,.5155,.5535,.5886/ data (absorp_table(i,9),i=1,30)/ & .0000,.0001,.0001,.0002,.0004,.0006,.0011,.0019,.0031,.0049, & .0075,.0110,.0157,.0221,.0309,.0430,.0591,.0798,.1056,.1366, & .1728,.2133,.2565,.3004,.3436,.3856,.4271,.4680,.5079,.5456/ data (absorp_table(i,8),i=1,30)/ & .0000,.0001,.0001,.0002,.0004,.0006,.0010,.0017,.0027,.0041, & .0060,.0086,.0122,.0170,.0236,.0323,.0437,.0585,.0772,.1005, & .1287,.1622,.2005,.2426,.2865,.3301,.3721,.4127,.4525,.4915/ data (absorp_table(i,7),i=1,30)/ & .0000,.0001,.0001,.0002,.0003,.0006,.0009,.0015,.0023,.0034, & .0049,.0070,.0099,.0138,.0189,.0256,.0342,.0454,.0597,.0778, & .1003,.1278,.1603,.1977,.2390,.2824,.3257,.3673,.4073,.4466/ data (absorp_table(i,6),i=1,30)/ & .0000,.0001,.0001,.0002,.0003,.0005,.0008,.0013,.0019,.0027, & .0039,.0056,.0079,.0109,.0149,.0200,.0265,.0349,.0457,.0595, & .0769,.0986,.1250,.1565,.1928,.2331,.2758,.3187,.3603,.4003/ data (absorp_table(i,5),i=1,30)/ & .0000,.0001,.0001,.0002,.0003,.0005,.0007,.0010,.0014,.0021, & .0029,.0042,.0058,.0079,.0107,.0143,.0189,.0247,.0321,.0415, & .0535,.0688,.0879,.1113,.1395,.1727,.2104,.2515,.2943,.3366/ data (absorp_table(i,4),i=1,30)/ & .0000,.0001,.0001,.0002,.0003,.0004,.0006,.0008,.0011,.0015, & .0021,.0029,.0039,.0052,.0069,.0091,.0118,.0152,.0194,.0248, & .0316,.0402,.0510,.0647,.0819,.1031,.1290,.1598,.1954,.2351/ data (absorp_table(i,3),i=1,30)/ & .0000,.0001,.0001,.0002,.0003,.0004,.0005,.0007,.0010,.0013, & .0018,.0024,.0032,.0042,.0054,.0069,.0087,.0109,.0137,.0171, & .0214,.0267,.0333,.0416,.0520,.0652,.0818,.1025,.1277,.1579/ data (absorp_table(i,2),i=1,30)/ & .0000,.0001,.0001,.0002,.0002,.0004,.0005,.0007,.0009,.0013, & .0017,.0022,.0029,.0037,.0046,.0058,.0071,.0086,.0105,.0127, & .0154,.0186,.0225,.0273,.0331,.0404,.0496,.0613,.0762,.0949/ data (absorp_table(i,1),i=1,30)/ & .0000,.0001,.0001,.0002,.0002,.0004,.0005,.0007,.0009,.0013, & .0017,.0022,.0028,.0036,.0044,.0054,.0066,.0079,.0094,.0111, & .0131,.0154,.0180,.0212,.0249,.0295,.0351,.0421,.0511,.0627/ data d_absorp_dt/ & 0.,.0006,.0017,.0036,.0061,.0092,.0158,0.,.0006,.0017,.0035,.0065, & .0114,.0230,0.,.0005,.0016,.0038,.0081,.0169,.0376,0.,.0005,.0017, & .0045,.0111,.0256,.0568,0.,.0004,.0021,.0065,.0178,.0421,.0864, & 0.,.0004,.0011,.0088,.0248,.0579,.1046,0.,.0004,.0032,.0112,.0316, & .0725,.1160,0.,.0004,.0036,.0141,.0392,.0880,.1256,0.,.0004,.0038, & .0187,.0508,.1045,.1345,0.,.0003,.0036,.0220,.0650,.1138,.1380,0., & .0002,.0028,.0234,.0859,.1180,.1410,0.,.0001,.0022,.0257,.0936, & .1183,.1444/ contains !####################################################################### subroutine mcm_lw_init(ix_in, jx_in, kx_in) integer, intent(in) :: ix_in, jx_in, kx_in ix = ix_in jx = jx_in kx = kx_in kp = kx + 1 km = kx -1 !------- write version number and namelist --------- if ( mpp_pe() == mpp_root_pe() ) then call write_version_number(version, tagname) endif module_is_initialized = .true. return end subroutine mcm_lw_init !####################################################################### subroutine mcm_lw_end module_is_initialized = .false. !--------------------------------------------------------------------- end subroutine mcm_lw_end !####################################################################### SUBROUTINE MCM_LW_RAD (KTOP,KBTM,NCLDS,EMCLD, & PRES,TEMP,RH2O,QO3,CAMT, & RRVCO2, HEATRA,GRNFLX,TOPFLX, phalf) ! This illustrates the intended sizes of arrays: ! INTEGER, INTENT(IN), DIMENSION(ix,jx,kp) :: KTOP,KBTM ! INTEGER, INTENT(IN), DIMENSION(ix,jx) :: NCLDS ! REAL, INTENT(IN), DIMENSION(ix,jx,kp) :: EMCLD ! REAL, INTENT(IN), DIMENSION(ix,jx,kp) :: PRES,TEMP ! REAL, INTENT(IN), DIMENSION(ix,jx,kx) :: RH2O,QO3 ! REAL, INTENT(IN), DIMENSION(ix,jx,kp) :: CAMT ! REAL, INTENT(IN) :: RRVCO2 ! REAL, INTENT(OUT), DIMENSION(ix,jx,kx) :: HEATRA ! REAL, INTENT(OUT), DIMENSION(ix,jx) :: GRNFLX,TOPFLX ! TK mod: ! REAL, INTENT(IN), DIMENSION(ix,jx,kp) :: phalf INTEGER, INTENT(IN), DIMENSION(:,:,:) :: KTOP,KBTM INTEGER, INTENT(IN), DIMENSION(:,:) :: NCLDS REAL, INTENT(IN), DIMENSION(:,:,:) :: EMCLD REAL, INTENT(IN), DIMENSION(:,:,:) :: PRES,TEMP REAL, INTENT(IN), DIMENSION(:,:,:) :: RH2O,QO3 REAL, INTENT(IN), DIMENSION(:,:,:) :: CAMT REAL, INTENT(IN) :: RRVCO2 REAL, INTENT(OUT), DIMENSION(:,:,:) :: HEATRA REAL, INTENT(OUT), DIMENSION(:,:) :: GRNFLX,TOPFLX ! TK mod: REAL, INTENT(IN), DIMENSION(:,:,:) :: phalf !----------------LOCAL ARRAY STORAGE------------------------------------ real, dimension(SIZE(PRES,1)) :: dummy REAL, DIMENSION(SIZE(PRES,1),SIZE(PRES,1),0:kp) :: cloud_cover REAL, DIMENSION(ix,jx,1:kx) :: sigma_level REAL, DIMENSION(ix,jx,0:kx) :: sigma_half_level REAL, DIMENSION(ix,jx,1:kx) :: sigma_thick REAL, DIMENSION(kx) :: o3_mixrat INTEGER :: i,j, n, ipr ! TK mods: integer :: iindex, jindex, klev, kpr real :: pi if(.not.module_is_initialized) then call error_mesg('MCM_LW_RAD','module is not initialized.',FATAL) endif ! ix = SIZE(PRES,1) ! jx = SIZE(PRES,2) lmax = ix * (kx + 1) ! Gather/calculate other needed parameters for ss longwave code: ! Convert grav to grav_accel in cgs units: grav_accel = grav * 100. pi = 4. * atan(1.) pi_alpha = 0.28 * pi t_freeze = tfreeze ! Multiply CO2 volume mixing ratio by 1.5194 to convert to mass mixing ! ratio... TK -- need to check accuracy of this... co2_mixrat = RRVCO2 * 1.5194 ! Compute sigma levels, which are needed by ss routine: ! This can be done on any x,y point in the domain. do klev = 1, kx do j=1,jx do i=1,ix sigma_level(i,j,klev) = pres(i,j,klev) / phalf(i,j,kp) enddo enddo enddo ! Compute half sigma levels, which are needed by ss routine: do klev = 1, kp do j=1,jx do i=1,ix sigma_half_level(i,j,klev-1) = phalf(i,j,klev) / phalf(i,j,kp) enddo enddo enddo ! Compute delta sigma of layers, which are needed by ss routine: do klev = 1, kx do j=1,jx do i=1,ix sigma_thick(i,j,klev) = sigma_half_level(i,j,klev) - & sigma_half_level(i,j,klev-1) enddo enddo enddo ! Create cloud_cover input field for ss longwave code: ! Expand_cloud routine creates a Manabe Climate Model cloud_cover ! field from the FMS arrays: nclds, ktop, kbtm, camt, and emcld. ! This is patterned after the expand_cloud subroutine in ! the clouds.f90 module. See also the supersource code ! in impt.f ! Accounts for cloud emmissivity... cloud_cover = 0.0 do jindex=1,size(nclds,2) do iindex=1,size(nclds,1) do n=2,nclds(iindex,jindex)+1 cloud_cover(iindex,jindex, & ktop(iindex,jindex,n):kbtm(iindex,jindex,n)) = & camt(iindex,jindex,n) * emcld(iindex,jindex,n) enddo enddo enddo ! TK Mod to mimic supersource. This is a temporary patch ! which needs to be cleaned up. A search is made for locations ! where the LW cloud emissivity has been set to 0.6. ! If that location is an isolated cloud with the layer above ! and below being cloud free, the cloud_cover value is left ! as 0.6, otherwise it is set to 1.0. This is the ! cirrus1L test. 7/09/01 do jindex=1,size(nclds,2) do iindex=1,size(nclds,1) ! For supersource runs, there should not be any cloud ! at the top model level (k=1), nor should the high cloud ! issue come into play the bottom level (k=kx). do klev = 2, kx-1 if (abs(cloud_cover(iindex,jindex,klev)-0.6) .lt. 0.001) then ! The cloud at level klev has been previously reset to 0.6 if ((cloud_cover(iindex,jindex,klev-1) .gt. 0.0) .or. & (cloud_cover(iindex,jindex,klev+1) .gt. 0.0)) then ! The cloud with 0.6 LW emiss has a neighboring cloud ! so reset the LW emiss (cloud_cover) to 1.0 as follows: cloud_cover(iindex,jindex,klev) = 1.0 end if end if enddo enddo enddo ! TK Can Artificially set cloud_cover to 1 at the test point: ! This method of doing this is now obsolete. See radiation_driver. ! print *, '***** TK Artificially set cloud_cover(49,27,14) = 1 ***' ! cloud_cover(49,27,14) = 1.0 DO j=1,jx ! Specify ozone for the latitude (zonal mean used): do klev=1,kx o3_mixrat(klev) = QO3(1,j,klev) end do ! TK Printout diagnostics: if ((j .eq. 27)) then ! if ((j .eq. 27) .or. (j .eq. 69)) then ! if ((j .eq. 10000) .or. (j .eq. 10000)) then ipr = 48 print * print *, 'TK chkpt 1 in mcm_lw_mod.f. Trap input to ' print *, ' lwcool for i,j = ', ipr, j print *, 'k cloud_cover PRESSURE' do kpr = 1, 14 write (6, 120) kpr, cloud_cover (ipr,j,kpr), & PRES(ipr,j,kpr) 120 format (i2, 2x, e10.3, 3x, e10.3) end do print * print *, 'Surface pressure = ', PRES(ipr,j,15) print * print *, 'k TEMP RH2O QO3 CAMT', & ' RRVCO2' do kpr = 1, 14 write (6, 125) kpr, TEMP(ipr,j,kpr),RH2O(ipr,j,kpr), & QO3(ipr,j,kpr), CAMT(ipr,j,kpr), RRVCO2 125 format (i2, 2x, f10.6, 2x, 2(e10.5,2x), 2(e8.3,2x)) end do print * print *, 'Surface temperature = ', TEMP(ipr,j,15) print * end if ! TK note that lw_down_sfc is just a throwaway variable in fms, ! so a dummy variable is used. Some extra input ! arguments are sent to lwcool. Send i=1 values of ! ozone since zonally symmetric. ! Use air temps in deg C as input parameter. ! Convert surface pressures from Pa to dynes/cm2 for ! input to the SS-derived routine... ! The extra parameter j is for debugging only... call lwcool ( cloud_cover(:,j,:), temp(:,j,1:kx)-t_freeze, & RH2O(:,j,1:kx), phalf(:,j,kp)*10., temp(:,j,kp), & dummy, HEATRA(:,j,:), TOPFLX(:,j), GRNFLX(:,j), j, & sigma_level(:,j,:), sigma_half_level(:,j,:), & sigma_thick(:,j,:), o3_mixrat) ! TK Printout diagnostics: if ((j .eq. 27) ) then ! if ((j .eq. 27) .or. (j .eq. 69)) then ! if ((j .eq. 10000) .or. (j .eq. 10000)) then ipr = 48 print * print *, 'TK chkpt 2b in mcm_lw_mod.f. Trap output from ' print *, ' lwcool for i,j = ', ipr, j ! Convert from ly/min to W/m2: print *, 'k HEATRA(K/d) GRNFLX TOPFLX (W/m2) ' do kpr = 1, 14 write (6, 130) kpr, HEATRA(ipr,j,kpr)*24*3600, & GRNFLX(ipr,j)*697., TOPFLX(ipr,j)*697. 130 format (i2, 2x, f10.6, 3x, f10.6, 3x, f10.6) end do print * end if ENDDO ! Convert top of atm and sfc fluxes into units FMS is ! expecting as output: (milliWatts/m2, K/day). ! TK mod 5/14/01: use identical conversion factor to SS: 697.496. ! Note that 4186/6 = 697.66667. GRNFLX = GRNFLX * 1.0E3 * 697.496 TOPFLX = TOPFLX * 1.0E3 * 697.496 ! GRNFLX = GRNFLX * 1.0E3 * 4186/6. ! TOPFLX = TOPFLX * 1.0E3 * 4186/6. HEATRA = HEATRA * 24 * 3600 END SUBROUTINE MCM_LW_RAD !####################################################################### subroutine lwcool(cloud_cover, tj, rj, psj, t_sfc, & lw_down_sfc, lw_cooling_rate, lw_up_toa, lw_net_sfc, jrow_flag, & sigma_level, sigma_half_level, sigma_thick, o3_mixrat) ! compute net longwave cooling rates ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! routine : lwcool ! called by : impt ! calls : lwtran ! purpose : compute net longwave cooling rates ! inputs: ! tj(i,k) temperature (deg C) at each grid point ! rj(i,k) mixing ratio of h2o (g/g) at each grid point ! psj(i) surface pressure (dynes/cm**2) at each grid point ! t_sfc(i) surface temperature (deg K) at each grid point ! cloud_cover(i,k) fractional cloud cover ! outputs: ! lw_up_toa(i) upward flux at top of atmosphere (ly/min) ! lw_net_sfc(i) net flux at ground (ly/min) ! lw_down_sfc(i) downward flux at ground (ly/min) ! lw_cooling_rate(i,k) cooling rate at each grid point (deg/s) ! procedure: ! the temperature at each grid point is converted from celsius to ! absolute (kelvin); from the temperatures planck blackbody ! radiances are computed for each grid point. the longwave ! transmission coefficients (recomputed every ldisk timesteps) are ! then combined with the blackbody radiances to determine clear ! fluxes (i.e., the flux which would occur in the absence of ! clouds) between all pairs of levels. dqx(i,k,l), for example, is ! the clear-path downward flux from vertical level l to level k. ! then, from the known distribution of clouds, the actual upward and ! downward fluxes at each level are computed. these are used to ! compute cooling rates as well as some diagnostic quantities. ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ real, parameter :: sigma = 5.673e-5, c1 = 3.740e-5, c2 = 1.4385 real, parameter :: cp1 = 0.24, cp2 = 4.1867e7 * cp1 real, intent(in) :: cloud_cover(ix,0:kp) real, intent(in) :: tj(ix,kx), rj(ix,kx), psj(ix) real, intent(in) :: t_sfc(ix) real, intent(in) :: sigma_level(ix,kx), sigma_half_level(ix,0:kx) real, intent(in) :: sigma_thick(ix,kx), o3_mixrat(kx) integer, intent(in) :: jrow_flag real, intent(out) :: lw_up_toa(ix) real, intent(out) :: lw_net_sfc(ix) real, intent(out) :: lw_cooling_rate(ix,kx) real, intent(out) :: lw_down_sfc(ix) !----------------LOCAL ARRAY STORAGE------------------------------------ real :: dqx (ix,0:kp,0:kp), uqx (ix,0:kp,0:kp) real :: dfxc(ix,0:kp), ufxc(ix,0:kp) real :: psigl(ix,2:kx) integer kloud (ix,0:kp) integer kldtop(ix,2:kx), kldmid(ix,2:kx), kldbot(ix,2:kx) integer :: numtop(2:kx), nummid(2:kx), numbot (2:kx) real :: ufxkb (ix), uqxkt (ix), ufxb(ix,2:kx) real :: dfxkt (ix), dqxkb (ix), dqxb(ix,2:kx) real :: ppfkb (ix), ppfkt (ix), ppfb(ix,2:kx) real :: cdxa (kx), cuxa(0:kx) real :: dfx(ix,0:kp), ufx(ix,0:kp) real temp_kelvin(ix,0:kp) real planck_func(ix,0:kp) real lw_trans_coeff(ix,0:kx,0:kx), lw_toa_correct(ix,0:kx) real lw_abs_quarter(ix,0:kx,2) real dpbb(ix,0:kx) integer :: k, l, kb, kt real :: cof ! TK: I took out these equivalences but needed to modify the ! code below so that dfx and ufx were initialized properly. ! Printout answers for test points reproduced. 5/03/01 ! equivalence (dfx, dqx(1,0,0)) ! equivalence (ufx, uqx(1,0,kp)) ! data cdxa / km * 0.5, 1.0 /, cuxa / -1.0, kx * -0.5 / do k=1,km cdxa(k) = 0.5 end do cdxa(kx) = 1.0 cuxa(0) = -1.0 do k=1,kx cuxa(k) = -0.5 end do ! ---------------------------------------------------------------------- ! compute clear (cloud-free) fluxes between all pairs of levels ! ---------------------------------------------------------------------- ! ****** compute pbb (planck blackbody radiance), and d(pbb)/dk do 5 k=1,kx do 5 i=1,ix temp_kelvin(i,k) = tj(i,k) + t_freeze 5 continue do 6 i=1,ix temp_kelvin(i, 0) = temp_kelvin(i,1) temp_kelvin(i,kp) = t_sfc(i) 6 continue do 7 k=0,kp do 7 i=1,ix planck_func(i,k) = sigma*temp_kelvin(i,k)**4 7 continue do 8 k=0,kx do 8 i=1,ix dpbb(i,k) = planck_func(i,k+1) - planck_func(i,k) 8 continue ! ---------------------------------------------------------------------- ! Obtain transmission functions ! ---------------------------------------------------------------------- call lwtran (planck_func, temp_kelvin, rj, psj, & lw_trans_coeff, lw_toa_correct, lw_abs_quarter, jrow_flag, & sigma_level, sigma_half_level, sigma_thick, o3_mixrat) ! ****** compute upward and downward fluxes at "flux" levels do 11 i=1,ix dqx(i,0,0) = 0.0 11 continue do 12 k=1,kx do 12 i=1,ix dqx(i,k,k) = cdxa(k) * dpbb(i,k) * lw_abs_quarter(i,k,2) + planck_func(i,k) 12 continue do 13 l=kx-1,0,-1 do 13 k=l+1,kx do 13 i=1,ix dqx(i,k,l) = dqx(i,k,l+1) - lw_trans_coeff(i,k,l)*dpbb(i,l) 13 continue do 14 k=1,kx do 14 i=1,ix dqx(i,k,0) = dqx(i,k,0) - lw_toa_correct(i,k) * planck_func(i,0) 14 continue do 15 k=0,kx do 15 i=1,ix uqx(i,k,k+1) = cuxa(k) * dpbb(i,k) * lw_abs_quarter(i,k,1) + planck_func(i,k+1) 15 continue do 16 l=2,kp do 16 k=0,l-2 do 16 i=1,ix uqx(i,k,l) = uqx(i,k,l-1) + lw_trans_coeff(i,k,l-1) * dpbb(i,l-1) 16 continue ! TK Add code to initialize ufx to take place of equivalence: do 17 k = 0, kp do 17 i = 1, ix ufx(i,k) = uqx(i,k,kp) 17 continue ! TK Add code to initialize dfx to take place of equivalence: do 18 k = 0, kp do 18 i = 1, ix dfx(i,k) = dqx(i,k,0) 18 continue ! ---------------------------------------------------------------------- ! compute actual fluxes from cloud distribution ! ---------------------------------------------------------------------- do 105 k=0,kp do 105 i=1,ix if (cloud_cover(i,k) .ne. 0) then kloud(i,k) = 1 else kloud(i,k) = 0 endif 105 continue do 110 k=2,kx do 110 i=1,ix kldtop(i,k) = kloud(i,k) * (1 - kloud(i,k-1)) kldmid(i,k) = kloud(i,k) * kloud(i,k+1) kldbot(i,k) = kloud(i,k) * (1 - kloud(i,k+1)) 110 continue do 20 kb=2,kx do 20 k=kb,kx do 20 i=1,ix if (kldbot(i,kb) .eq. 1) & dfx(i,k) = dfx(i,k) * (1.0-cloud_cover(i,kb)) & + dqx(i,k,kb) * cloud_cover(i,kb) 20 continue do 30 kt=kx,2,-1 do 30 k=0,kt-1 do 30 i=1,ix if (kldtop(i,kt) .eq. 1) & ufx(i,k) = ufx(i,k) * (1.0-cloud_cover(i,kt)) & + uqx(i,k,kt) * cloud_cover(i,kt) 30 continue do 35 k=2,kx numtop(k) = 0 numbot(k) = 0 nummid(k) = 0 35 continue do 36 k=2,kx do 36 i=1,ix numtop(k) = numtop(k) + kldtop(i,k) numbot(k) = numbot(k) + kldbot(i,k) nummid(k) = nummid(k) + kldmid(i,k) 36 continue do 40 k=kx,2,-1 do 115 i=1,ix if (kldbot(i,k) .eq. 1) then ufxkb(i) = ufx(i,k) dqxkb(i) = dqx(i,k,k) ppfkb(i) = sigma_half_level(i,k) endif 115 continue if (nummid(k) .eq. 0) goto 40 do 120 i=1,ix if (kldmid(i,k) .eq. 1) then ufxb(i,k) = ufxkb(i) dqxb(i,k) = dqxkb(i) ppfb(i,k) = ppfkb(i) endif 120 continue 40 continue do 50 k=2,kx if (nummid(k) .eq. 0) goto 50 do 125 i=1,ix if (kldtop(i,k) .eq. 1) then dfxkt(i) = dfx(i,k-1) uqxkt(i) = uqx(i,k-1,k) ppfkt(i) = sigma_half_level(i,k-1) endif 125 continue do 130 i=1,ix if (kldmid(i,k) .eq. 1) then psigl(i,k) = (sigma_half_level(i,k) - ppfkt(i)) / & (ppfb(i,k) - ppfkt(i)) dfxc (i,k) = dfxkt(i) + psigl(i,k)*(dqxb(i,k) - dfxkt(i)) ufxc (i,k) = uqxkt(i) + psigl(i,k)*(ufxb(i,k) - uqxkt(i)) endif 130 continue 50 continue do 135 k=2,kx do 135 i=1,ix if (kldmid(i,k) .eq. 1) then dfx(i,k) = dfx(i,k)*(1.0-cloud_cover(i,k)) + dfxc(i,k)*cloud_cover(i,k) ufx(i,k) = ufx(i,k)*(1.0-cloud_cover(i,k)) + ufxc(i,k)*cloud_cover(i,k) endif 135 continue ! ---------------------------------------------------------------------- ! compute net cooling rates from fluxes ! ---------------------------------------------------------------------- ! 1.43306e-6 converts ergs / (cm**2*sec) to ly/min 80 continue do 85 i=1,ix lw_up_toa(i) = ufx(i,0) * 1.43306e-06 lw_net_sfc(i) = (ufx(i,kx) - dfx(i,kx)) * 1.43306e-06 lw_down_sfc(i) = -dfx(i,kx) * 1.43306e-6 85 continue ! 1.1574e-5 = 1.0 / 86400.0 cof = grav_accel * 86400.0 / cp2 cof = cof * 1.1574e-5 do 90 k=1,kx do 90 i=1,ix lw_cooling_rate(i,k) = cof/(sigma_thick(i,k)*psj(i)) 90 continue do 95 k=1,kx do 95 i=1,ix lw_cooling_rate(i,k) = lw_cooling_rate(i,k) & * (ufx(i,k) - ufx(i,k-1) + dfx(i,k-1) - dfx(i,k)) 95 continue return end subroutine lwcool function expx(x) real, intent(in) :: x real :: expx ! ---------------------------------------------------------------------- ! statement function ! ---------------------------------------------------------------------- ! the approximation to exp(x) which follows is appropriate only given ! the following assumptions: ! 1) since the correction coefficients (h2o_corr_a & h2o_corr_b) ! are known only to about four digits of accuracy, 6-7 digits ! of accuracy for the approximation is ample. ! 2) the input range for x is [-6.735,2.289]. this depends on: ! a) the valid temperature range, now [-173.16,102] (celsius) ! b) the values of the coefficients h2o_corr_a & h2o_corr_b ! real, parameter :: c0e = 9.999999810120017e-1 real, parameter :: c1e = 9.999999599140620e-1/16.**1 real, parameter :: c2e = 5.000056658124412e-1/16.**2 real, parameter :: c3e = 1.666837874572320e-1/16.**3 real, parameter :: c4e = 4.146341773620255e-2/16.**4 real, parameter :: c5e = 7.279860860195404e-3/16.**5 expx = ((( (((((c5e*x+c4e)*x+c3e)*x+c2e)*x+c1e)*x+c0e) & **2)**2)**2)**2 end function expx subroutine lwtran (planck_func, temp_kelvin, rj, psj, & lw_trans_coeff, lw_toa_correct, lw_abs_quarter, jrow_flag, & sigma_level, sigma_half_level, sigma_thick, o3_mixrat) ! compute longwave atmospheric transmission coefficients ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! routine : lwtran ! called by : lwcool ! calls : lwco2 ! purpose : compute longwave atmospheric transmission coefficients ! inputs: ! rj(i,k) mixing ratio of h2o at each grid point (g/g) ! psj(i,k) surface presure at each ground point (dynes/cm**2) ! temp_kelvin(i,k) temperature (deg K) ! planck_func(i,k) blackbody function (total radiance; ergs/(cm**2)*s) ! outputs: ! lw_trans_coeff longwave transmission coefficient (dimensionless) ! lw_toa_correct correction to above for top level (dimensionless) ! lw_abs_quarter quarter-level absorption (dimensionless) ! procedure: ! the main purpose of this routine is to compute the band-integrated ! transmission coefficients transm(i,k,l) for radiation passing from ! flux level l to flux level k. in addition, it computes two ! corrections. the first, transz(i,k), is for radiation from the top ! level (split off from transm to improve code efficiency, i think). ! the second, absrbx(i,k,kq), is a half-layer (flux level to data- ! input level, or vice versa) absorption correction; kq=1 is the ! correction for upward flux and kq=2 is for downward flux. ! this routine splits the longwave radiation spectrum into nb bands, ! ranging from about 20 to 2000 cm-1. transmission coefficients are ! computed for each absorber contributing to absorption in that ! band. the single-absorber transmission coefficients are ! multiplied together to produce the total single-band transmission, ! and then subtracted from 1 to produce an absorption. the single- ! band absorptions are weighted by the relative fraction of ! blackbody radiation in that band (for a particular temperature), ! and summed to produce the total band-integrated absorption ! coefficient, which is converted back to a transmission. ! transmissions for h2o are computed using a goody random model. ! the function form and coefficients come from a paper by rogers and ! walshaw (rogers, c.d. and c.d. walshaw, the computation of infra- ! red cooling in planetary atmospheres, quart. j. royal meteorol. ! soc., v.92, pp.67-92, 1966). ! the water vapor continuum is now computed by using the ! formulation given in roberts(first fit). reference: ! roberts,r.e.,et al,applied optics,v. 15,pp 2085-2090,1976. ! the co2 transmission table was generated by a program obtained ! privately from s.r. drayson, university of michigan. this program ! assumed a homogeneous path between flux levels; a version which ! did not assume homogeneous paths was later published by drayson ! (drayson, s.r., atmospheric transmission in the co2 bands between ! 12 microns and 18 microns, appl. opt., v.5, pp.385-391, 1973). ! the o3 transmission computation, changed from a table lookup in ! the original longwave radiation code, also uses a goody random ! model. the functional form and coefficients were published by ! rogers (rogers, c.d., some extensions and applications of the new ! random model for molecular band transmission, q.j.r.m.s, v.94, ! pp.99-102, 1968). ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ real, intent(in) :: planck_func(ix,0:kp) real, intent(in) :: temp_kelvin(ix,0:kp) real, intent(in) :: rj(ix,kx) real, intent(in) :: psj(ix) real, intent(in) :: sigma_level(ix,kx), sigma_half_level(ix,0:kx) real, intent(in) :: sigma_thick(ix,kx), o3_mixrat(kx) integer, intent(in) :: jrow_flag real, intent(out) :: lw_trans_coeff(ix,0:kx,0:kx) real, intent(out) :: lw_toa_correct(ix,0:kx) real, intent(out) :: lw_abs_quarter(ix,0:kx,2) real, parameter :: sigma = 5.673e-5, c1 = 3.740e-5, c2 = 1.4385 real, parameter :: po = 1013.25e3, opo = 1.0 / po, o2po = 0.5 * opo ! parameter ( opo = 9.8692271e-7, o2po = 4.9346136e-7 ) real, parameter :: ak1 = 223.0, ak2=10.0 ! parameter ( alpha = 0.28 ) ! parameter ( pi = 3.1415926535898 ) ! parameter ( pialph = pi * alpha ) real, parameter :: ak1x4 = 4.0 * ak1 real, parameter :: ak2x4 = 4.0 * ak2 ! parameter ( pialph = 0.8796459, ak1x4 = 892.0, ak2x4 = 40.0 ) real, parameter :: azint = 1.66 ! parameter ( g = 980.6, azint = 1.66, g1 = azint / g ) ! parameter ( g1h = 0.5 * g1 ) ! parameter ( g1 = 1.6938774e-3, g1h = 0.5 * g1 ) real, parameter :: c3 = 1.0 / 120.0, a1 = 75.48, a2 = 23.44 ! parameter ( c3 = 8.333333e-3, a1 = 75.48, a2 = 23.44 ) real, parameter :: c3a1 = c3 * a1, c3a2 = c3 * a2 real, parameter :: c0e = 9.999999810120017e-1 real, parameter :: c1e = 9.999999599140620e-1/16.**1 real, parameter :: c2e = 5.000056658124412e-1/16.**2 real, parameter :: c3e = 1.666837874572320e-1/16.**3 real, parameter :: c4e = 4.146341773620255e-2/16.**4 real, parameter :: c5e = 7.279860860195404e-3/16.**5 real :: bw(ix,0:kx,nb), bwz(ix,nb), expzs(ix,nb) real :: trans(ix,0:kx,0:kx), dbbdt(ix,0:kx) real :: vcube(nb), c(nb), v(nb), dv(nb) real :: pdi (ix,0:kp), pfl (ix,0:kx), dpdi (ix, kx) real :: rh2o (ix,0:kp), ro3 (ix,0:kp), qdi (ix, kx) real :: duh2o (ix, kx), duco2 (ix, kx), duo3 (ix, kx) real :: duqh2o(ix, kx), duqfac(ix, kx), dufac(ix, kx) real :: fbignl(ix, kx), bignel(ix,0:kx), bign (ix, kx) real :: tdpfac(ix, kx), tfac (ix, kx), teff (ix, kx) real :: sfac (ix,0:kx), pfac (ix,0:kx), peff (ix,0:kx) real :: tcof (ix,0:kx), fac1 (ix,0:kx), fac2 (ix,0:kx) real :: t260 (ix,0:kx), tadj (ix,0:kx,2) real :: ulog (ix,0:kx), plog (ix,0:kx) ! quarter-layer variables real :: txdegk(ix,0:kx), qx (ix,0:kx), dpfac(ix,0:kx) real :: duxh2o(ix,0:kx), duxco2(ix,0:kx), duxo3(ix,0:kx) real :: g1, g1h, bc, bandc1, bandc2, cco2, asodsq, diag, & factr1, factr2, diago3 integer :: lb, k, l, len, kq !c TK: Removed unnecessary equivalences: ! equivalence (ulog, fac1 ), (plog, fac2 ), (teff, peff ) ! equivalence (txdegk,dpdi ), (qx, qdi ) ! equivalence (duxh2o,duh2o ), (duxco2,duco2), (duxo3,duo3 ) ! equivalence (tadj, duqfac), (tadj(1,0,2),dufac,tdpfac,dpfac) ! constants for h2o bignell continuum absorption ! these are the roberts continuum coefficients data c / 4*0.0, 149.566, 38.103, 11.312, 7.9193, 6.0741, 4.9179, 9*0.0 / ! longwave frequency bands and band widths (inverse cm) data v / 20.0, 100.0, 220.0, 340.0, 470.0, & 670.0, 850.0, 930.0, 1020.0, 1140.0, & 1275.0, 1400.0, 1500.0, 1600.0, 1700.0, & 1800.0, 1900.0, 2000.0, 2125.0 / data dv / 40.0, 120.0, 120.0, 120.0, 140.0, & 260.0, 100.0, 60.0, 120.0, 120.0, & 150.0, 100.0, 100.0, 100.0, 100.0, & 100.0, 100.0, 100.0, 150.0 / ! ---------------------------------------------------------------------- ! description of variables and constants used in lwtran ! ---------------------------------------------------------------------- ! variables centered on data input levels: ! sigma_level 1-d "sigma" (pressure/po) ! sigma_thick d(sigma)/dk ! (not to be confused with planck bb sigma) ! pdi, dpdi, qdi pressure, d(pdi)/dk, and pdi/po ! rh2o, co2_mixrat, ro3 mixing ratios for the various absorbers ! duh2o, duco2, duo3 absorber amounts from top of atmosphere ! (gm/cm**2) (also called unadjusted optical depth) ! temp_kelvin temperature (kelvin) ! variables centered on flux levels (see last paragraph, below): ! sigma_half_level, pfl 1-d "sigma", 2-d pressure ! assorted additional variables and constants: ! po, opo, o2po average pressure at sea level, 1/po, 2/po ! planck_func planck blackbody function ! dbbdt d(planck_func)/dt ! sigma, c1, c2 constants in planck blackbody equation ! v, dv frequencies and band widths (1/cm) ! (wavelength range is 5-500 microns) ! bw, bwz weights for freq. bands derived from planck_func ! alpha, c3, a1, a2, empirical constants used in computing ! ak1, ak2 ozone absorption ! h2o_line_width, ! h2o_line_strength empirical band constants for water lines ! h2o_line_width coefficient for h2o line width ! h2o_line_strength coefficient for h2o line strength ! sfac, pfac, tfac intensity, pressure, temperature factors ! peff, teff effective pressure, temperature ! ulog, plog logs of co2 amount, effective pressure ! tadj temperature adjustments for h2o lines ! fbignl, bignel, bign factors associated with bignell continuum ! tcof transmission coefficient, single absorber ! trans transmission matrix, single frequency band ! output variables: ! transm, transz band-integrated transmission coefficients ! absrbx quarter-layer correction to above (absrp.) ! quarter layer variables (x usually indicates quarter-layer): ! txdegk, qx quarter-layer temperature, "sigma" ! duxh2o, duxco2, duxo3 quarter-layer absorber amounts ! miscellaneous constants: ! azint azimuthal integral factor for diffuse rad. ! g1, g1h convert dp*mix-ratio to absorber amount ! miscellaneous temporaries: ! expzs, vcube, bc, bandc1, bandc2 ! t260, fac1, fac2, dufac, dpfac, duqfac, duqh2o, cco2 ! vertical levels in the model are numbered from the top of the ! atmosphere to the bottom; thus vertical level 1 is in the ! stratosphere and level kx is just above the ground. in the ! radiation code, additional "pseudo-levels" are tacked onto the top ! (level 0) and the bottom (level kp). these levels make it easier ! to handle incoming solar radiation and longwave radiation emanating ! from the ground. ! an additional complication to the level numbering system for the ! radiation code is the use of "flux" levels and "data-input" levels. ! the simplest numbering scheme for radiation calculations focuses ! on even and odd levels: the grid variables (pressure, temperature, ! etc.) reside at even levels, and flux is computed across them, i.e., ! from odd level to odd level. this a pedagogically sound approach, ! but it poses problems for computers, particularly vector types. ! thus this code numbers what would be odd and even levels separately. ! "data-input" levels are normal model vertical levels; "flux" levels ! are halfway between them. flux level k is one-half level closer to ! the ground than data-input level k. in the diagram below, data-input ! levels are shown as "++++++" and flux levels are shown as "------". ! description level index ! incoming solar radiation 0 0 +++++++ ! 1/2 0 ------- ! top model level 1 1 +++++++ ! . ! . ! . ! bottom model level kx kx +++++++ ! kx+1/2 kx ------- ! ground level kp kp +++++++ ! ---------------------------------------------------------------------- ! Compute various constants involving gravity. ! ---------------------------------------------------------------------- g1 = azint / grav_accel g1h = 0.5 * g1 ! ---------------------------------------------------------------------- ! compute band weights ! ---------------------------------------------------------------------- do 50 lb=1,nb vcube(lb) = v(lb) * v(lb) ** 2 50 continue do 52 lb=1,nb do 52 i=1,ix expzs(i,lb) = exp(c2 / temp_kelvin(i,0) * v(lb)) - 1.0 52 continue do 54 lb=1,nb bc = c1 * vcube(lb) * dv(lb) do 54 i=1,ix bwz(i,lb) = bc / (expzs(i,lb) * planck_func(i,0)) 54 continue do 56 k=0,kx do 56 i=1,ix dbbdt(i,k) = 4.0 * sigma * temp_kelvin(i,k)**3 56 continue do 58 lb=1,nb bandc1 = c1 * vcube(lb) bandc2 = c2 * v(lb) bc = bandc1 * bandc2 * dv(lb) do 58 k=0,kx do 58 i=1,ix bw(i,k,lb) = exp(bandc2 / temp_kelvin(i,k)) bw(i,k,lb) = bc * bw(i,k,lb) & / ((bw(i,k,lb) - 1.0) * temp_kelvin(i,k))**2 bw(i,k,lb) = bw(i,k,lb) / dbbdt(i,k) 58 continue ! ---------------------------------------------------------------------- ! compute field arrays needed for trans calculation ! (note that temp_kelvin has already been computed and stored in flux) ! ---------------------------------------------------------------------- do 70 i=1,ix pdi(i, 0) = 0.0 pfl(i, 0) = 0.0 pdi(i,kp) = psj(i) 70 continue do 72 k=1,kx do 72 i=1,ix pdi (i,k) = psj(i) * sigma_level(i,k) pfl (i,k) = psj(i) * sigma_half_level(i,k) dpdi(i,k) = psj(i) * sigma_thick(i,k) ro3 (i,k) = o3_mixrat(k) 72 continue do 74 k=1,kx do 74 i=1,ix rh2o (i,k) = rj(i,k) qdi (i,k) = (pfl(i,k) + pfl(i,k-1)) * o2po t260 (i,k) = temp_kelvin(i,k) - 260.0 fbignl(i,k) = exp(1800.0/temp_kelvin(i,k)-6.081081081) ! the kernel is 1800.0*(1.0/temp-1.0/296.0) 74 continue do 76 i=1,ix rh2o(i, 0) = 0.0 ro3 (i, 0) = 0.0 rh2o(i,kp) = rh2o(i,kx) ro3 (i,kp) = ro3 (i,kx) 76 continue do 77 k=1,kp do 77 i=1,ix rh2o(i,k) = max(rh2o(i,k),3.0e-6) 77 continue cco2 = co2_mixrat * g1 do 78 k=1,kx do 78 i=1,ix duh2o (i,k) = dpdi (i,k) * rh2o (i,k) * g1 duco2 (i,k) = dpdi (i,k) * cco2 duo3 (i,k) = dpdi (i,k) * ro3 (i,k) * g1 duqh2o(i,k) = duh2o (i,k) * qdi (i,k) bignel(i,k) = (fbignl(i,k)*duqh2o(i,k)*rh2o(i,k)) / (0.622+rh2o(i,k)) 78 continue ! ---------------------------------------------------------------------- ! compute transm and transz ! transm(i,k,l) is the transmission coefficient for longwave radiation ! passing from vertical level l to level k ! transz is a correction for radiation from the top of the atmosphere ! ---------------------------------------------------------------------- do 80 l=0,kx do 80 k=0,kx do 80 i=1,ix trans (i,k,l) = 0.0 lw_trans_coeff(i,k,l) = 0.0 80 continue do 82 k=0,kx do 82 i=1,ix lw_toa_correct(i,k) = 0.0 82 continue ! ***** begin loop on longwave frequency bands ***** do 2890 lb= 1, nb ! ****** h2o contribution ! ****** compute temperature adjustment to line strength and line width do 1300 l=1,2 do 1300 k=1,kx do 1300 i=1,ix tadj(i,k,l) = (h2o_corr_a(lb,l) + h2o_corr_b(lb,l) * t260(i,k)) * t260(i,k) tadj(i,k,l) = expx(tadj(i,k,l)) 1300 continue do 1310 k=0,kx do 1310 i=1,ix sfac(i,k) = 0.0 pfac(i,k) = 0.0 1310 continue do 1320 k=1,kx do 1320 i=1,ix bign(i,k) = 0.0 1320 continue ! compute h2o contribution to lower triangular half of trans matrix ! note that a non-physical factor of sod has been folded into duqfac ! to cancel the factor that arises when sfac is squared in 1550 loop asodsq = h2o_line_strength(lb) * h2o_line_width(lb) if (h2o_line_strength(lb) .eq. 0.0) asodsq = h2o_line_width(lb) do 1400 k=1,kx do 1400 i=1,ix dufac (i,k) = h2o_line_strength(lb) * duh2o (i,k) * tadj(i,k,2) duqfac(i,k) = asodsq * duqh2o(i,k) * tadj(i,k,1) 1400 continue do 1550 l=kx-1,0,-1 do 1500 k=l+1,kx do 1500 i=1,ix sfac(i,k) = sfac(i,k) + dufac (i,l+1) pfac(i,k) = pfac(i,k) + duqfac(i,l+1) bign(i,k) = bign(i,k) + bignel(i,l+1) 1500 continue do 1550 k=l+1,kx do 1550 i=1,ix trans(i,k,l) = exp(-c(lb) * bign(i,k) - sfac(i,k) / sqrt(1.0 + sfac(i,k)**2 / pfac(i,k))) 1550 continue ! ****** compute h2o contribution to diagonal elements of trans matrix diag = exp(-c(lb) - h2o_line_strength(lb)/sqrt(1.0 + h2o_line_strength(lb)/h2o_line_width(lb))) ! ****** co2 contribution if (lb .ne. 6) goto 2500 do 2005 k=0,kx do 2005 i=1,ix sfac(i,k) = 0.0 pfac(i,k) = 0.0 2005 continue do 2007 k=1,kx do 2007 i=1,ix tfac(i,k) = 0.0 2007 continue ! **** compute co2 contribution to lower triangular half of trans matrix do 2010 k=1,kx do 2010 i=1,ix duqfac(i,k) = duco2(i,k) * qdi (i,k) tdpfac(i,k) = temp_kelvin(i,k) * dpdi(i,k) 2010 continue do 2050 l=kx-1,0,-1 do 2020 k=l+1,kx do 2020 i=1,ix sfac(i,k) = sfac(i,k) + duco2 (i,l+1) pfac(i,k) = pfac(i,k) + duqfac(i,l+1) tfac(i,k) = tfac(i,k) + tdpfac(i,l+1) teff(i,k) = tfac(i,k) / (pfl(i,k) - pfl(i,l)) 2020 continue do 2030 k=l+1,kx do 2030 i=1,ix ulog(i,k) = alog10(sfac(i,k)) plog(i,k) = alog10(pfac(i,k)) - ulog(i,k) ulog(i,k) = ulog(i,k) + 2.70668 2030 continue ! ------ ORIGINAL CODE ------------------------ ! if (l .eq. kx-1) then ! i = 0 ! ulog(i,kx) = 2.70668 ! plog(i,kx) = 0.0 ! teff(i,kx) = 0.0 ! len = ix + 1 ! else ! i = 1 ! len = ix * (kx - l) ! endif ! ****** call lwco2 to compute tcof values ! call lwco2(ulog(i,l+1), plog(i,l+1), teff(i,l+1), len, tcof(i,l+1)) ! ****** compute co2 contribution to diagonal elements of trans matrix ! if (l .eq. kx-1) diag = tcof(i,kx) * diag ! ---------END OF ORIGINAL CODE ----------------------- ! ---------BEGIN TK MODIFICATIONS if (l .eq. kx-1) then i = ix ulog(i,kx-1) = 2.70668 plog(i,kx-1) = 0.0 teff(i,kx-1) = 0.0 len = ix + 1 ! ****** call lwco2 to compute tcof values call lwco2(ulog(i,kx-1), plog(i,kx-1), teff(i,kx-1), & len, tcof(i,kx-1)) ! ****** compute co2 contribution to diagonal elements of trans matrix diag = tcof(i,kx-1) * diag else i = 1 len = ix * (kx - l) ! ****** call lwco2 to compute tcof values call lwco2(ulog(i,l+1), plog(i,l+1), teff(i,l+1), & len, tcof(i,l+1)) endif ! -----------END OF TK MODIFICATIONS ------------------- ! ****** accumulate co2 tcof into trans (lower triangular half) do 2040 k=l+1,kx do 2040 i=1,ix trans(i,k,l) = tcof(i,k) * trans(i,k,l) 2040 continue 2050 continue ! ****** ozone contribution 2500 if (lb .ne. 9) go to 2855 do 2555 k=0,kx do 2555 i=1,ix sfac(i,k) = 0.0 pfac(i,k) = 0.0 2555 continue ! sfac = 0.0 ! pfac = 0.0 ! **** compute o3 contribution to lower triangular half of trans matrix do 2600 k=1,kx do 2600 i=1,ix duqfac(i,k) = pi_alpha * duo3(i,k) * qdi(i,k) 2600 continue do 2700 l=kx-1,0,-1 do 2650 k=l+1,kx do 2650 i=1,ix sfac(i,k) = sfac(i,k) + duo3 (i,l+1) pfac(i,k) = pfac(i,k) + duqfac(i,l+1) 2650 continue do 2700 k=l+1,kx do 2700 i=1,ix peff (i,k ) = pfac(i,k) / sfac(i,k) fac1 (i,k ) = sfac(i,k) / peff(i,k) fac2 (i,k ) = 1.0 - exp(-5.0 * peff(i,k) * (sqrt(1.0 + ak2x4 * fac1(i,k)) - 1.0)) fac1 (i,k ) = 1.0 - exp(-5.0 * peff(i,k) * (sqrt(1.0 + ak1x4 * fac1(i,k)) - 1.0)) tcof (i,k ) = 1.0 - (c3a1 * fac1(i,k) + c3a2 * fac2(i,k)) trans(i,k,l) = tcof(i,k) * trans(i,k,l) 2700 continue ! ****** compute o3 contribution to diagonal elements of trans matrix factr1 = 1.0 - exp(-5.0 * pi_alpha * (sqrt(1.0 + ak1x4 / pi_alpha) - 1.0)) factr2 = 1.0 - exp(-5.0 * pi_alpha * (sqrt(1.0 + ak2x4 / pi_alpha) - 1.0)) diago3 = 1.0 - (c3a1 * factr1 + c3a2 * factr2) diag = diago3 * diag ! ****** apply band weights to trans; accumulate ! ****** as absorp. in transm & transz ! ****** first convert trans and diag to absorption 2855 continue do 2857 l=0,kx do 2857 k=0,kx do 2857 i=1,ix trans(i,k,l) = 1.0 - trans(i,k,l) 2857 continue diag = 1.0 - diag do 2860 l=0,kx-1 do 2860 k=l+1,kx do 2860 i=1,ix lw_trans_coeff(i,k,l) = lw_trans_coeff(i,k,l) + trans(i,k,l) * bw(i,l,lb) 2860 continue do 2865 l=1,kx do 2865 k=0,l-1 do 2865 i=1,ix lw_trans_coeff(i,k,l) = lw_trans_coeff(i,k,l) + trans(i,l,k) * bw(i,l,lb) 2865 continue do 2870 l=0,kx do 2870 i=1,ix lw_trans_coeff(i,l,l) = lw_trans_coeff(i,l,l) + diag * bw(i,l,lb) 2870 continue do 2875 l=1,kx do 2875 i=1,ix lw_toa_correct(i,l) = lw_toa_correct(i,l) + trans(i,l,0) * bwz(i,lb) 2875 continue do 2880 i=1,ix lw_toa_correct(i,0) = lw_toa_correct(i,0) + diag * bwz(i,lb) 2880 continue 2890 continue ! ***** end loop on longwave frequency bands ***** ! convert transm & transz from absorption to transmission coefficients do 2897 l=0,kx do 2897 k=0,kx do 2897 i=1,ix lw_trans_coeff(i,k,l) = 1.0 - lw_trans_coeff(i,k,l) 2897 continue do 2898 k=0,kx do 2898 i=1,ix lw_toa_correct(i,k) = 1.0 - lw_toa_correct(i,k) 2898 continue ! ---------------------------------------------------------------------- ! compute absrbx ! absrbx is a correction for the finite thickness of the vertical ! levels, and represents within-layer absorption of radiation ! ---------------------------------------------------------------------- ! absrbx = 0.0 do 2899 l=1,2 do 2899 k=0,kx do 2899 i=1,ix lw_abs_quarter(i,k,l) = 0.0 2899 continue ! kq=1 is the upward-flux correction; kq=2 is for downward flux do 5000 kq=1,2 ! compute coefficients needed for quarter-layer correction if (kq.eq.2) goto 3050 do 3010 k=1,kx-1 do 3010 i=1,ix dpfac(i,k) = (pdi(i,k+1) - pfl(i,k)) * g1h 3010 continue do 3020 k=1,kx-1 do 3020 i=1,ix txdegk(i,k) = 0.375 * temp_kelvin(i,k) + 0.625 * temp_kelvin(i,k+1) qx (i,k) = (0.75 * pfl(i,k) + 0.25 * pdi(i,k+1)) * opo duxh2o(i,k) = (0.375 * rh2o(i,k) + 0.625 * rh2o(i,k+1)) * dpfac(i,k) duxco2(i,k) = co2_mixrat * dpfac(i,k) duxo3 (i,k) = (0.375 * ro3(i,k) + 0.625 * ro3(i,k+1)) * dpfac(i,k) bignel(i,k) = (0.375 * rh2o(i,k) + 0.625 * rh2o(i,k+1))**2 & * dpfac(i,k) * fbignl(i,k+1) / & (0.622 + 0.375*rh2o(i,k) + 0.625*rh2o(i,k+1)) 3020 continue do 3030 i=1,ix txdegk(i, 0) = 0.75 * temp_kelvin(i,0) + 0.25 * temp_kelvin(i,1) txdegk(i,kx) = temp_kelvin(i,kp) qx (i, 0) = (0.75 * pdi(i,0) + 0.25 * pdi(i,1)) * opo qx (i,kx) = psj(i) * opo duxh2o(i, 0) = (0.25 * rh2o(i,0) + 0.75 * rh2o(i,1)) * (pdi(i,1) - pdi(i,0)) * g1h duxh2o(i,kx) = 0.0 duxco2(i, 0) = co2_mixrat * (pdi(i,1) - pdi(i,0)) * g1h duxco2(i,kx) = 0.0 duxo3 (i, 0) = (0.25 * ro3(i,0) + 0.75 * ro3(i,1)) * (pdi(i,1) - pdi(i,0)) * g1h duxo3 (i,kx) = 0.0 bignel(i, 0) = (0.75 * rh2o(i,0) + 0.25 * rh2o(i,1))**2 & * (pdi(i,1)-pdi(i,0)) * g1h * fbignl(i,1) / & (0.622 + 0.75*rh2o(i,0) + 0.25*rh2o(i,1)) bignel(i,kx) = 0.0 3030 continue do 3040 k=0,kx do 3040 i=1,ix t260(i,k) = txdegk(i,k) - 260.0 3040 continue go to 3100 3050 do 3060 k=1,kx-1 do 3060 i=1,ix dpfac(i,k) = (pfl(i,k) - pdi(i,k)) * g1h 3060 continue do 3070 k=1,kx-1 do 3070 i=1,ix txdegk(i,k) = 0.625 * temp_kelvin(i,k) + 0.375 * temp_kelvin(i,k+1) qx (i,k) = (0.75 * pfl(i,k) + 0.25 * pdi(i,k)) * opo duxh2o(i,k) = (0.625 * rh2o(i,k) + 0.375 * rh2o(i,k+1)) * dpfac(i,k) duxco2(i,k) = co2_mixrat * dpfac(i,k) duxo3 (i,k) = (0.625 * ro3(i,k) + 0.375 * ro3(i,k+1)) * dpfac(i,k) bignel(i,k) = (0.625*rh2o(i,k) + 0.375*rh2o(i,k+1))**2 & * dpfac(i,k) * fbignl(i,k) / & (0.622 + 0.625*rh2o(i,k) + 0.375*rh2o(i,k+1)) 3070 continue do 3080 i=1,ix txdegk(i, 0) = temp_kelvin(i,0) txdegk(i,kx) = 0.75 * temp_kelvin(i,kp) + 0.25 * temp_kelvin(i,kx) qx (i, 0) = 1.0 qx (i,kx) = (0.75 * pdi(i,kp) + 0.25 * pdi(i,kx)) * opo duxh2o(i, 0) = 0.0 duxh2o(i,kx) = (0.75 * rh2o(i,kp) + 0.25 * rh2o(i,kx)) * (pdi(i,kp) - pdi(i,kx)) * g1h duxco2(i, 0) = 0.0 duxco2(i,kx) = co2_mixrat * (pdi(i,kp) - pdi(i,kx)) * g1h duxo3 (i, 0) = 0.0 duxo3 (i,kx) = (0.75 * ro3(i,kp) + 0.25 * ro3(i,kx)) * (pdi(i,kp) - pdi(i,kx)) * g1h bignel(i, 0) = 0. bignel(i,kx) = (0.75 * rh2o(i,kp) + 0.25 * rh2o(i,kx))**2 & * (pdi(i,kp)-pdi(i,kx)) * g1h * fbignl(i,kx) / & (0.622 + 0.75*rh2o(i,kp) + 0.25*rh2o(i,kx)) 3080 continue do 3090 k=0,kx do 3090 i=1,ix t260(i,k) = txdegk(i,k) - 260.0 3090 continue ! ***** begin loop on longwave frequency bands ***** 3100 do 4000 lb=1,nb ! ****** quarter-layer h2o absorption ! ****** compute temperature adjustment to line strength and line width do 3150 k=0,kx do 3150 i=1,ix tadj(i,k,1) = (h2o_corr_a(lb,1) + h2o_corr_b(lb,1) * t260(i,k)) * t260(i,k) tadj(i,k,2) = (h2o_corr_a(lb,2) + h2o_corr_b(lb,2) * t260(i,k)) * t260(i,k) tadj(i,k,1) = tadj(i,k,1) - tadj(i,k,2) 3150 continue do 3200 l=1,2 do 3200 k=0,kx do 3200 i=1,ix tadj(i,k,l) = expx(tadj(i,k,l)) 3200 continue ! ****** compute h2o contribution to quarter-layer absorption matrix do 3300 k=0,kx do 3300 i=1,ix sfac (i,k ) = h2o_line_strength(lb) * duxh2o(i,k) * tadj(i,k,2) pfac (i,k ) = h2o_line_width(lb) * qx(i,k) * tadj(i,k,1) trans(i,k,kq) = exp(-c(lb) * bignel(i,k) -sfac(i,k)/sqrt(1.0+sfac(i,k)/pfac(i,k))) 3300 continue ! ****** quarter-layer co2 absorption if (lb .ne. 6) go to 3700 ! ****** compute co2 contribution to quarter-layer absorption matrix len = ix * kp do 3410 k=0,kx do 3410 i=1,ix sfac(i,k) = duxco2(i,k) 3410 continue do 3420 k=0,kx do 3420 i=1,ix duxco2(i,k) = max(duxco2(i,k),1.0e-6) 3420 continue do 3500 k=0,kx do 3500 i=1,ix ulog(i,k) = alog10(duxco2(i,k)) + 2.70668 plog(i,k) = alog10(qx(i,k)) 3500 continue ! ****** call lwco2 to compute tcof values call lwco2(ulog(1,0), plog(1,0), txdegk(1,0), len, tcof(1,0)) ! ****** accumulate co2 ccp into trans do 3600 k=0,kx do 3600 i=1,ix tcof(i,k) = tcof(i,k) * trans(i,k,kq) 3600 continue do 3610 k=0,kx do 3610 i=1,ix if(sfac(i,k) .gt. 0.0) then trans(i,k,kq) = tcof(i,k) endif 3610 continue ! ****** quarter-layer o3 absorption 3700 if (lb .ne. 9) go to 3900 do 3800 k=0,kx do 3800 i=1,ix peff (i,k ) = pi_alpha * qx(i,k) fac2 (i,k ) = duxo3(i,k) / peff(i,k) fac1 (i,k ) = 1.0 - exp(-5.0 * peff(i,k) * (sqrt(1.0 + ak1x4 * fac2(i,k)) - 1.0)) fac2 (i,k ) = 1.0 - exp(-5.0 * peff(i,k) * (sqrt(1.0 + ak2x4 * fac2(i,k)) - 1.0)) tcof (i,k ) = 1.0 - (c3a1 * fac1(i,k) + c3a2 *fac2(i,k)) trans(i,k,kq) = tcof(i,k) * trans(i,k,kq) 3800 continue ! ****** total quarter-layer absorption 3900 do 4000 k=0,kx do 4000 i=1,ix lw_abs_quarter(i,k,kq) = lw_abs_quarter(i,k,kq) + (1.0 - trans(i,k,kq)) * bw(i,k,lb) 4000 continue 5000 continue ! ***** end of principal quarter-layer loop ***** return end subroutine lwtran subroutine lwco2 (ulog,plog,teff,ln,tco2) ! table lookup routine for co2 transmission ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! routine: lwco2 ! called by: lwtran ! calls: none ! purpose: compute longwave co2 transmission ! inputs: ! ulog log of co2 amount ! plog log of effective pressure ! teff effective temperature ! ln length of call-argument arrays ! output: ! tco2 co2 transmission coefficient ! side effects: ! none ! (i.e., no input variables or common blocks modified, etc.) ! procedure: ! a 2d table lookup of co2 absorption is performed. the axes are ! log of the co2 amount and log of the effective pressure. an ! effective temperature of 300k is assumed for this first step. ! a correction for temperatures in the range 200k-300k is made by ! multiplying a temperature correction factor by the difference ! between the effective temperature and 300k, and adding the ! result to the previously calculated absorption. the temperature ! correction factor is found by 2d table lookup, using the same ! axes as were used to obtain the 300k absorption. no temperature ! correction is applied if the temperature exceeds 300k; if the ! temperature falls below 200k, the 200k correction is applied. ! if the co2 concentration is too low, no temperature correction ! is made. ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! TK Mod: set up in parent routine: ! parameter ( lmax = ix * (kx + 1) ) integer, intent(in) :: ln real, intent(in) :: ulog(ln) real, intent(in) :: plog(ln) real, intent(in) :: teff(ln) real, intent(out) :: tco2(ln) real :: ahi (lmax), alo (lmax), cct (lmax), cdt (lmax), fulog (lmax) integer :: nplog (lmax), iulog (lmax), index (lmax) real :: plnum (lmax), plden (lmax), plfac (lmax), tcor (lmax) real :: alogoo(lmax), alogpo(lmax), alogop(lmax), alogpp(lmax) real :: pltab ( 12) integer :: iu, l, n !c TK: Removed unnecessary equivalences: !c TK: equivalence (alogoo, plnum), (alogpo, plden, alo, tcor) !c TK: equivalence (alogop, iulog), (alogpp, index, ahi, cdt ) data pltab / & -4.50000, -4.00000, -3.40000, -2.88081, -2.27875, -1.88081, & -1.57978, -1.27875, -0.88081, -0.48287, 0. , 0.29528 / ! ****** check for out-of-range ulog input values iu = 0 do 25 l=1,ln if ( ulog(l).lt.-4.25 ) iu=iu+1 25 continue if ( iu.ne.0 ) go to 1200 ! ****** compute pressure range for tables, ! ****** and pressure interpolation factor do 50 l=1,ln nplog(l) = 1 50 continue do 100 n=2,11 do 100 l=1,ln if ( plog(l).gt.pltab(n) ) then nplog(l) = n plnum(l) = plog(l) - pltab(n) plden(l) = pltab(n+1) - pltab(n) endif 100 continue do 200 l=1,ln plfac(l) = plnum(l) / plden(l) 200 continue ! ****** compute co2 absorption for an effective temperature of 300k do 500 l=1,ln fulog(l) = 4 * ulog(l) + 18 iulog(l) = int(fulog(l)) ! TK diagnostic printout: ! print *, 'TK chkpt 1 in lwco2: l, ulog(l), fulog(l), ', & ! 'int(fulog(l)), iulog(l) = ', & ! l, ulog(l), fulog(l), int(fulog(l)), iulog(l) if ( iulog(l).gt.29 ) iulog(l) = 29 fulog(l) = fulog(l) - real(iulog(l)) 500 continue do 600 l = 1, ln index(l) = iulog(l) + 30 * (nplog(l) - 1) 600 continue do 620 l=1,ln alogoo(l) = absorp_table(index(l),1) 620 continue do 640 l=1,ln index (l) = index(l) + 1 alogpo(l) = absorp_table(index(l),1) 640 continue do 660 l=1,ln index (l) = index(l) + (30 - 1) alogop(l) = absorp_table(index(l),1) 660 continue do 680 l=1,ln index (l) = index(l) + 1 alogpp(l) = absorp_table(index(l),1) 680 continue do 700 l = 1, ln ahi(l) = alogop(l) + fulog(l) * (alogpp(l) - alogop(l)) alo(l) = alogoo(l) + fulog(l) * (alogpo(l) - alogoo(l)) cct(l) = alo (l) + plfac(l) * (ahi (l) - alo (l)) 700 continue ! ****** compute temperature correction coefficient for co2 absorption do 750 l=1,ln fulog(l) = ulog(l) + 4 iulog(l) = int(fulog(l)) if ( iulog(l).gt.6 ) iulog(l) = 6 fulog(l) = fulog(l) - real(iulog(l)) 750 continue do 800 l = 1, ln index(l) = iulog(l) + 7 * (nplog(l) - 1) 800 continue do 820 l=1,ln alogoo(l) = d_absorp_dt(index(l),1) 820 continue do 840 l=1,ln index (l) = index(l) + 1 alogpo(l) = d_absorp_dt(index(l),1) 840 continue do 860 l=1,ln index (l) = index(l) + (7 - 1) alogop(l) = d_absorp_dt(index(l),1) 860 continue do 880 l=1,ln index (l) = index(l) + 1 alogpp(l) = d_absorp_dt(index(l),1) 880 continue do 900 i = 1, ln ahi(i) = alogop(i) + fulog(i) * (alogpp(i) - alogop(i)) alo(i) = alogoo(i) + fulog(i) * (alogpo(i) - alogoo(i)) cdt(i) = alo (i) + plfac(i) * (ahi (i) - alo (i)) 900 continue ! ****** restrict effective temperature input values to proper range do 950 l=1,ln tcor(l) = teff(l) if ( tcor(l).gt.320.0 ) tcor(l) = 320 if ( tcor(l).lt.180.0 ) tcor(l) = 180 950 continue ! ****** apply temperature correction to co2 ! ****** absorption if ulog is in range do 1000 l=1,ln cdt (l) = cct(l) + 0.01 * (cdt(l) * (tcor(l) - 300)) tco2(l) = 1 - 1.3077 * cdt(l) 1000 continue do 1100 l=1,ln if ( ulog(l).lt.-3 ) tco2(l) = 1 - 1.3077 * cct(l) 1100 continue return ! ****** fatal error: ulog out-of-range 1200 write(6,1201) ulog call Error_Mesg ( 'lwco2 in MCM_LW_MOD','ulog out of range.', FATAL) 1201 format('1ulog > 3. or < -4.25; execution terminated' / ' the ulog vector is:' / (8e16.8)) end subroutine lwco2 end module mcm_lw_mod