c $Header: /fs/cgd/csm/models/CVS.REPOS/atm/ccm_crm_src/crm/Attic/crm.F,v 1.1.2.14 1999/09/07 22:33:02 zender Exp $
c CCM3 Column Radiation Model
c Jeffrey Kiehl, Bruce Briegleb, and Charlie Zender
c July 1998
c All routines except crm(), getdat(), netcdf(), radctl(), radcsw() and
c four dummy routines (readric(), writeric(), outfld(), and radozn())
c are #included straight from CCM code. Thus CRM code is identical to CCM code
c but does not invoke the Land Surface Model for land albedo computation.
c Routines radctl() and radcsw() are slightly modified to support the CRM.
c The purpose of the (non-dummy) routines is:
c crm(): The main() routine. This routine takes the place of tphysbc()
c in the CCM3 calling tree. It places the required calls to the cloud
c routines cldefr() and cldems() directly, rather than invoking cldint().
c getdat(): Reads the stdin file to parse the user-specified column
c profile. Overwrites the co2vmr previsouly set by radini() with the
c user specified value. Computes o3vmr (instead of in radozn()).
c Sets calday variable in comtim.h used in zenith() and radinp().
c radctl(): Main radiation driving routine, same as ccm3 version except:
c receives additional input variable o3vmr, and passes out
c additional diagnostic radiation quantities and eccf usually local to radctl().
c does cgs-->mks conversion for all fluxes.
c netcdf(): Routine to output a netCDF file
c named crm.nc containing most of the data.
c radcsw(): Altered to compute SRB quantities and pass them into crmsrb.h
#include <misc.h>
#include <params.h>
c The following block allows the entire CRM to compile with a simple command
c Using this command also checks that all needed files are in the CRM
#ifdef SINGLE_SOURCE_FILE
#include <aermix.F>
#include <albocean.F>
#include <blkdat.F>
#include <cldefr.F>
#include <cldems.F>
#include <endrun.F>
#include <freemem.F>
#include <getmem.F>
#include <fmrgrid.F>
#include <orb.F>
#include <radabs.F>
#include <radclr.F>
#include <radclw.F>
#include <radcsw.F>
#include <radctl.F>
#include <radded.F>
#include <radems.F>
#include <radini.F>
#include <radinp.F>
#include <radoz2.F>
#include <radtpl.F>
#include <resetr.F>
#include <torgrid.F>
#include <trcab.F>
#include <trcabn.F>
#include <trcems.F>
#include <trcmix.F>
#include <trcplk.F>
#include <trcpth.F>
#include <zenith.F>
#include <netcdf.F>
#ifndef CRAY
#include <intmax.F>
#include <isrchfgt.F>
#include <isrchfle.F>
#include <wheneq.F>
#include <whenfgt.F>
#include <whenflt.F>
#include <whenne.F>
#endif /* CRAY */
#endif /* not SINGLE_SOURCE_FILE */
program crm,8
#include <implicit.h>
c Parameters
#include <prgrid.h>
c Commons
#include <comtim.h> /* calday */
#include <comvmr.h> /* co2vmr, n2ovmr, ch4vmr, cfc11, cfc12 */
#include <comsol.h> /* scon, tauvis, eccen, obliq, mvelp, iyear_AD */
#include <comctl.h> /* anncyc,iradsw,iradlw,iradae */
#ifdef CRM_SRB
#include <crmsrb.h>
#endif /* not CRM_SRB */
c Arguments
c Fields specified by the user in getdat()
real clon(plon) ! Centered longitude (radians)
real clat ! Current centered latitude (radians)
real cld(plond,plevp) ! fractional cloud cover
real clwp(plond,plev) ! cloud liquid water path
real coslat ! cosine latitude
c NB: o3mmr and o3vmr should be dimensioned (plond,plevr) if a different
c size radiation grid is used. Clashes between prgrid.h and ptrrgrid.h
c (they both define plngbuf) prevent us from dimensioning anything by
c plevr in this top level crm() routine.
real o3mmr(plond,plev) ! Ozone mass mixing ratio
real o3vmr(plond,plev) ! Ozone volume mixing ratio
real aldif(plond) ! Albedo: longwave, diffuse
real aldir(plond) ! Albedo: longwave, direct
real asdif(plond) ! Albedo: shortwave, diffuse
real asdir(plond) ! Albedo: shortwave, direct
real oro(plond) ! Land/ocean/sea ice flag
real pilnm1(plond,plevp) ! natural log of pintm1
real pintm1(plond,plevp) ! model interface pressures
real pmidm1(plond,plev) ! model level pressures
real pmlnm1(plond,plev) ! natural log of pmidm1
real ps(plond) ! surface pressure
real qm1(plond,plev) ! model level specific humidity
real snowh(plond) ! snow depth (liquid water equivalent)
real tg(plond) ! surface (skin) temperature
real tm1(plond,plev) ! model level temperatures
real ts(plond) ! surface air temperature
c Fields computed from user input
real coszrs(plond) ! cosine solar zenith angle
real eccf ! earth/sun distance factor
real effcld(plond,plevp) ! effective cloud=cld*emis
real emis(plond,plev) ! cloud emissivity
real fice(plond,plev) ! fractional amount of ice
real loctim(plond) ! local time of solar computation
real lwup(plond) ! Longwave up flux at surface
real rei(plond,plev) ! ice particle size
real rel(plond,plev) ! liquid effective drop size (microns)
real srfrad(plond) ! srf radiative heat flux
c Output longwave arguments from radctl()
real flwds(plond) ! Surface down longwave flux
real qrl(plond,plev) ! Longwave cooling rate
c Output shortwave arguments from radctl()
real fsns(plond) ! Surface absorbed solar flux
real qrs(plond,plev) ! Solar heating rate
real soll(plond) ! Downward solar rad onto surface (lw direct)
real solld(plond) ! Downward solar rad onto surface (lw diffuse)
real sols(plond) ! Downward solar rad onto surface (sw direct)
real solsd(plond) ! Downward solar rad onto surface (sw diffuse)
c Additional CRM diagnostic output from radctl()
real flns(plond) ! srf longwave cooling (up-dwn) flux
real flnsc(plond) ! clr sky lw flx at srf (up-dwn)
real flnt(plond) ! net outgoing lw flx at model top
real flntc(plond) ! clr sky lw flx at model top
real fsnsc(plond) ! clr sky surface abs solar flux
real fsnt(plond) ! total column absorbed solar flux
real fsntc(plond) ! clr sky total column abs solar flux
real solin(plond) ! solar incident flux
real fsnirt(plond) ! [W m-2] Near-IR flux absorbed at TOA
real fsnirtsq(plond) ! [W m-2] Near-IR flux absorbed at TOA>= 0.7 microns
real fsnrtc(plond) ! [W m-2] Clear sky near-IR flux absorbed at TOA
c Local workspace: These variables are not saved
real hbuf ! history buffer
real pie ! 3.14159...
integer i ! longitude index
integer k ! level index
integer lat ! latitude row index
c Fundamental constants needed by radini()
real cpair ! heat capacity dry air at constant prs (J/kg/K)
real epsilo ! ratio mean mol weight h2o to dry air
real gravit ! gravitational acceleration (m/s**2)
real stebol ! Stefan-Boltzmann constant (W/m**2/K**4)
c Externals
external blkdat
c Main Code
#ifdef SUN
c 19990907: IEEE code deprecated by default
c CCM IEEE handling code from control/ccm3.F
c#include <f77_floatingpoint.h>
c integer iexcept
c integer ieee_handler
c iexcept=ieee_handler('set','common',SIGFPE_ABORT)
c if (iexcept.ne.0) write(6,*) 'IEEE trapping not supported here'
#endif /* not SUN */
write(6,'(a)') 'Begin CCM3 Column Radiation Model'
c Set latitude index to 1 (only used in calling radctl())
c CCM: dynamics/eul/scan1bc()
lat=1
c Set control/comtim.h: nstep,dosw,dolw,doabsems
c CCM: dynamics/eul/stepon()
nstep=0
c CCM: dynamics/advnce()
dosw=.true.
c CCM: dynamics/advnce()
dolw=.true.
doabsems=.true.
c Set control/comctl.h: anncyc,iradsw,iradlw,iradae
c CCM: control/data()
anncyc=.true.
c CCM: control/data()
iradsw=1
c CCM: control/data()
iradlw=1
c CCM: control/data()
iradae=1
c Set physics/comcon.h: gravit, cpair, epsilo, stebol
c NB: physics/comcon.h is never actually used in the CRM
c These SI parameters are used in physics/inti() to call physics/radini() to set CGS parameters of same name in physics/crdcon.h: gravit, cpair, epsilo, stebol
c These constants must be set before computing lwup or calling radini()
c CCM: dynamics/eul/initcom()
gravit = 9.80616
cpair = 1004.64
epsilo = 0.622
stebol = 5.67e-8
c Besides setting the variables passed as subroutines arguments, getdat() sets
c control/comtim.h: calday (needed in zenith() and radinp())
c physics/comsol.h: scon, tauvis, iyear_AD, obliqr, lambm0 and mvelpp
c physics/comvmr.h: co2vmr, n2ovmr, ch4vmr, f11vmr, f12vmr
c control/crdcon.h: pie
call getdat(
$ aldif,
$ aldir,
$ asdif,
$ asdir,
$ clat,
$ cld,
$ clon,
$ clwp,
$ coslat,
$ loctim,
$ o3mmr,
$ o3vmr,
$ oro,
$ pilnm1,
$ pintm1,
$ pmidm1,
$ pmlnm1,
$ ps,
$ qm1,
$ snowh,
$ tg,
$ tm1,
$ ts)
c Given these four physical constants in MKS, radini() defines CGS equivalents,
c and sets many radiation parameters stored in common blocks.
c radini() must be called after getdat(), because the CO2 mass mixing ratio
c set in radini() depends on co2vmr set by the user in getdat().
c CCM: physics/inti()
call radini(gravit,cpair,epsilo,stebol)
c Get coszrs: needed as input to albocean() and radctl()
c CCM: dom/initext()
call zenith (calday ,clat ,clon ,coszrs )
c Adjust the user-specified albedo for ocean/sea-ice points only
c Land points always use the user-specified albedos
c CCM: dom/initext()
call albocean(oro ,snowh ,coszrs ,asdir ,aldir ,
$ asdif ,aldif )
c Cloud particle size and fraction of ice
c CCM: physics/cldint()
call cldefr(oro, tm1, rel, rei, fice, ps, pmidm1)
c Cloud emissivity
c CCM: physics/cldint()
call cldems(clwp, fice, rei, emis)
c The radiation scheme recomputes the skin temperature in physics/radtpl()
c This skin temperature is based on the LW up flux at the surface
c Thus, surface LW up flux is computed here, then passed down to radtpl() where,
c hopefully, it is inverted to obtain the same skin temperature we started with.
c CCM: control/initext() (Ocean and sea-ice only)
c CCM: ? (Land)
do i=1,plon
lwup(i)=stebol*tg(i)**4 ! W m-2 LW up is skin temperature to fourth power
end do
c Effective cloud cover
c CCM: physics/cldint()
do k=1,plev
do i=1,plon
effcld(i,k) = cld(i,k)*emis(i,k)
end do
end do
c Cloud cover at surface interface always zero (for safety's sake)
c CCM: physics/cldint()
do i=1,plon
effcld(i,plevp) = 0.0
cld(i,plevp) = 0.0
end do
c Main radiation driving routine
c NB: All fluxes returned from radctl() have already been converted to MKS
c CCM: physics/tphysbc()
call radctl(hbuf ,clat ,coslat ,lat ,lwup ,
$ pmidm1 ,pintm1 ,pmlnm1 ,pilnm1 ,tm1 ,
$ qm1 ,cld ,effcld ,clwp ,coszrs ,
$ asdir ,asdif ,aldir ,aldif ,fsns ,
$ qrs ,qrl ,flwds ,rel ,rei ,
$ fice ,sols ,soll ,solsd ,solld
c++csz
$ ,fsnt,fsntc,fsnsc,flnt,flns,flntc,flnsc,solin, ! Output
$ fsnirt,fsnirtsq,fsnrtc, ! Output
$ eccf, ! Output (computed in radctl()/radinp()/orbdecl())
$ o3vmr) ! Input
c--csz
c CCM: physics/tphysbc()
do i=1,plon
srfrad(i)=fsns(i)+flwds(i)
end do
#ifdef CRM_SRB
c Diagnostic column abundances
do i=1,plon
mpc_H2O(i)=0.0 ! kg m-2
mpc_O3(i)=0.0 ! kg m-2
end do
do i=1,plon
do k=1,plev
pdelm1(i,k)=pintm1(i,k+1)-pintm1(i,k) ! Pa
mpc_H2O(i)=mpc_H2O(i)+pdelm1(i,k)*qm1(i,k)/gravit ! kg m-2
mpc_O3(i)=mpc_O3(i)+pdelm1(i,k)*o3mmr(i,k)/gravit ! kg m-2
end do ! end loop over lev
mpc_O3_DU(i)=mpc_O3(i)*gas_cst_O3*tpt_STP/prs_STP ! m
mpc_O3_DU(i)=mpc_O3_DU(i)*1.0e5 ! m --> millicm = 1.0e-5 m
end do ! end loop over lon
#endif /* not CRM_SRB */
c Write out final results
write (6,'(a)') 'CCM3 CRM Results:'
write (6,'(a)') 'Conventions:'
write (6,*) 'Shortwave fluxes are positive downward'
write (6,*) 'Longwave fluxes are positive upward'
write (6,*) 'Net Radiative fluxes are positive downward (into the system)'
write (6,*) 'Fluxes defined to be zero are not reported (e.g., LW down flx TOA)'
write (6,'(a)') 'Abbreviations, Acronyms and Definitions:'
write (6,*) 'LW = Longwave'
write (6,*) 'LWCF = Longwave Cloud Forcing'
write (6,*) 'NCF = Net Cloud Forcing = SWCF+LWCF'
write (6,*) 'NIR = Near Infrared (0.7 < lambda < 5.0 microns)'
write (6,*) 'N7 = NOAA7 satellite NIR instrument weighted flux'
write (6,*) 'NRF = Net Radiative Flux: sum of SW and LW fluxes'
write (6,*) 'SW = Shortwave'
write (6,*) 'SWCF = Shortwave Cloud Forcing'
write (6,*) 'TOA = Top of Atmosphere'
write (6,*) 'Vis = Visible (0.2 < lambda < 0.7 microns)'
write (6,*) 'atm = Atmosphere'
write (6,*) 'clr = Clear sky (diagnostic computation with no clouds)'
write (6,*) 'ctr = Center'
write (6,*) 'dff = Diffuse flux'
write (6,*) 'drc = Direct flux'
write (6,*) 'dwn = Downwelling'
write (6,*) 'frc = Fraction'
write (6,*) 'lqd = Liquid'
write (6,*) 'mpc = Mass path column'
write (6,*) 'net = Net flux = downwelling minus upwelling flux'
write (6,*) 'spc = Spectral'
write (6,*) 'sfc = Surface level'
write (6,*) 'vmr = Volume mixing ratio'
write (6,*) 'wvl = Wavelength'
write (6,*) 'um = Microns'
write (6,*) 'up = Upwelling'
write (6,*) ''
write (6,'(a)') 'Sun-Earth Geometry:'
pie = 4.0*atan(1.0) ! NB: This is not the pie in physics/crdcon.h
do i=1,1
write (6,*) 'Year AD = ',iyear_AD
write (6,*) 'Day of year (Greenwich) = ',calday
write (6,*) 'Local solar hour = ',loctim(i)
write (6,*) 'Latitude = ',180.0*clat/pie,' degrees'
write (6,*) 'Longitude = ',180.0*clon(i)/pie,' degrees'
write (6,*) 'Solar zenith angle = ',180.0*acos(coszrs(i))/pie, ' degrees'
write (6,*) 'Cosine solar zenith angle = ',coszrs(i)
write (6,*) 'Earth-sun distance = ',eccf,' AU'
write (6,*) 'Solar constant = ',scon/1000.0,' W m-2' ! CGS --> SI
write (6,*) ''
write (6,'(a)') 'Shortwave (SW) results ( < 5.0 um):'
if (solin(i).le.0.0) then
write (6,*) 'SW albedo = Ill-defined'
write (6,*) 'SW albedo (clr) = Ill-defined'
else
write (6,*) 'SW albedo = ',(solin(i)-fsnt(i))/solin(i)
write (6,*) 'SW albedo (clr) = ',(solin(i)-fsntc(i))/solin(i)
endif
write (6,*) 'SW down flux TOA = ',solin(i),' W m-2'
write (6,*) 'SW up flux TOA = ',solin(i)-fsnt(i),' W m-2'
write (6,*) 'SW up flux TOA (clr) = ',solin(i)-fsntc(i),' W m-2'
write (6,*) 'SW net flux TOA = ',fsnt(i),' W m-2'
write (6,*) 'SW net flux TOA (clr) = ',fsntc(i),' W m-2'
write (6,*) 'SW flux abs atm = ',fsnt(i)-fsns(i),' W m-2'
write (6,*) 'SW flux abs atm (clr) = ',fsntc(i)-fsnsc(i),' W m-2'
write (6,*) 'SW down flux sfc = ',sols(i)+soll(i)+solsd(i)+solld(i),' W m-2'
c write (6,*) 'SW down flux sfc (clr) = ',,' W m-2'
write (6,*) 'SW up flux sfc = ',sols(i)+soll(i)+solsd(i)+solld(i)-fsns(i),' W m-2'
c write (6,*) 'SW up flux sfc (clr) = ',,' W m-2'
write (6,*) 'SW net flux sfc = ',fsns(i),' W m-2'
write (6,*) 'SW net flux sfc (clr) = ',fsnsc(i),' W m-2'
write (6,*) 'SW cloud forcing TOA = ',fsnt(i)-fsntc(i),' W m-2'
write (6,*) 'SW cloud forcing sfc = ',fsns(i)-fsnsc(i),' W m-2'
if (fsnt(i)-fsntc(i).eq.0.0) then
write (6,*) 'SWCF(sfc)/SWCF(TOA) = Ill-defined'
else
write (6,*) 'SWCF(sfc)/SWCF(TOA) = ',(fsns(i)-fsnsc(i))/(fsnt(i)-fsntc(i))
endif ! endif
write (6,*) ''
write (6,'(a)') 'Longwave (LW) results ( > 5.0 um):'
write (6,*) 'LW up flux TOA = ',flnt(i),' W m-2'
write (6,*) 'LW up flux TOA (clr) = ',flntc(i),' W m-2'
write (6,*) 'LW up flux sfc = ',flns(i)+flwds(i),' W m-2'
write (6,*) 'LW down flux sfc = ',flwds(i),' W m-2'
c write (6,*) 'LW down flux sfc (clr) = ',flwds(i),' W m-2'
write (6,*) 'LW net flux sfc = ',flns(i),' W m-2'
write (6,*) 'LW net flux sfc (clr) = ',flnsc(i),' W m-2'
write (6,*) 'LW cloud forcing TOA = ',flntc(i)-flnt(i),' W m-2'
write (6,*) 'LW cloud forcing sfc = ',flnsc(i)-flns(i),' W m-2'
write (6,*) ''
write (6,'(a)') 'Net Radiative Flux results (NRF=SW+LW):'
write (6,*) 'NRF up flux TOA = ',solin(i)-fsnt(i)+flnt(i),' W m-2'
write (6,*) 'NRF down flux TOA = ',solin(i),' W m-2'
c write (6,*) 'NRF up flux TOA (clr) = ',flntc(i),' W m-2'
write (6,*) 'NRF net flux TOA = ',fsnt(i)-flnt(i),' W m-2'
write (6,*) 'NRF net flux TOA (clr) = ',fsntc(i)-flntc(i),' W m-2'
write (6,*) 'NRF up flux sfc = ',sols(i)+soll(i)+solsd(i)+solld(i)-fsns(i)+flns(i)+flwds(i),' W m-2'
write (6,*) 'NRF down flux sfc = ',sols(i)+soll(i)+solsd(i)+solld(i)+flwds(i),' W m-2'
c write (6,*) 'NRF down flux sfc (clr) = ',flwds(i),' W m-2'
write (6,*) 'NRF net flux sfc = ',fsns(i)-flns(i),' W m-2'
write (6,*) 'NRF net flux sfc (clr) = ',fsnsc(i)-flnsc(i),' W m-2'
write (6,*) 'NRF cloud forcing TOA = ',fsnt(i)-fsntc(i)+flntc(i)-flnt(i),' W m-2'
write (6,*) 'NRF cloud forcing sfc = ',fsns(i)-fsnsc(i)+flnsc(i)-flns(i),' W m-2'
write (6,*) ''
#ifdef CRM_SRB
write (6,'(a)') 'Specified atmospheric constituents:'
write (6,*) 'Visible AOD = ',tauvis
write (6,*) 'H2O mpc = ',mpc_H2O(i),' kg m-2'
write (6,*) 'O3 mpc = ',mpc_O3(i),' kg m-2'
write (6,*) 'O3 mpc = ',mpc_O3_DU(i),' Dobson'
write (6,*) 'CO2 vmr = ',co2vmr
write (6,*) 'N2O vmr = ',n2ovmr
write (6,*) 'CH4 vmr = ',ch4vmr
write (6,*) 'F11 vmr = ',f11vmr
write (6,*) 'F12 vmr = ',f12vmr
write (6,*) ''
write (6,'(a)') 'Column extinction optical depths:'
write (6,'(1x,a,f6.4,a,f6.4,a)')
$ 'Visible band = ',wvl_min(bnd_idx_vsb)*1.0e6,'--',wvl_max(bnd_idx_vsb)*1.0e6,' um'
write (6,*) 'Tau total = ',odxc_ttl(i,bnd_idx_vsb)
write (6,*) 'Tau Ray = ',odxc_Ray(i,bnd_idx_vsb)
write (6,*) 'Tau aer = ',odxc_aer(i,bnd_idx_vsb)
write (6,*) 'Tau lqd = ',odxc_lqd(i,bnd_idx_vsb)
write (6,*) 'Tau ice = ',odxc_ice(i,bnd_idx_vsb)
write (6,*) 'Tau O3 = ',odxc_O3(i,bnd_idx_vsb)
write (6,*) 'Tau H2O = ',odxc_H2O(i,bnd_idx_vsb)
write (6,*) 'Tau O2 = ',odxc_O2(i,bnd_idx_vsb)
write (6,*) 'Tau CO2 = ',odxc_CO2(i,bnd_idx_vsb)
write (6,*) ''
write (6,'(a)') 'Visible spectral fluxes:'
write (6,'(1x,a,f6.4,a,f6.4,a)')
$ 'Visible band = ',wvl_min(bnd_idx_vsb)*1.0e6,'--',wvl_max(bnd_idx_vsb)*1.0e6,' um'
write (6,*) 'Down spc flux TOA = ',flx_bnd_dwn_TOA(i,bnd_idx_vsb)*1.0e-6/wvl_dlt(bnd_idx_vsb), ' W m-2 um-1'
write (6,*) 'Up spc flux TOA = ',flx_bnd_up_TOA(i,bnd_idx_vsb)*1.0e-6/wvl_dlt(bnd_idx_vsb), ' W m-2 um-1'
write (6,*) 'Down spc flux sfc = ',flx_bnd_dwn_sfc(i,bnd_idx_vsb)*1.0e-6/wvl_dlt(bnd_idx_vsb), ' W m-2 um-1'
write (6,*) 'Down spc flux dff sfc = ',flx_bnd_dwn_dff_sfc(i,bnd_idx_vsb)*1.0e-6/wvl_dlt(bnd_idx_vsb), ' W m-2 um-1'
write (6,*) 'Down spc flux drc sfc = ',flx_bnd_dwn_drc_sfc(i,bnd_idx_vsb)*1.0e-6/wvl_dlt(bnd_idx_vsb), ' W m-2 um-1'
write (6,*) 'Up spc flux sfc = ',flx_bnd_up_sfc(i,bnd_idx_vsb)*1.0e-6/wvl_dlt(bnd_idx_vsb), ' W m-2 um-1'
write (6,*) ''
c Output TOA Radiation Budget (TRB)
write (6,'(a)') 'Solar TOA radiation budget components:'
if (solin(i).le.0.0) then
write (6,*) 'SW alb TOA = Ill-defined'
write (6,*) 'Vis alb TOA = Ill-defined'
write (6,*) 'NIR alb TOA = Ill-defined'
write (6,*) 'alb(NIR)/alb(SW) TOA = Ill-defined'
write (6,*) 'alb(NIR)/alb(Vis) TOA = Ill-defined'
else
write (6,*) 'SW alb TOA = ',alb_SW_TOA(i)
write (6,*) 'Vis alb TOA = ',alb_vsb_TOA(i)
write (6,*) 'NIR alb TOA = ',alb_NIR_TOA(i)
write (6,*) 'alb(NIR)/alb(SW) TOA = ',alb_NIR_SW_TOA(i)
write (6,*) 'alb(NIR)/alb(Vis) TOA = ',alb_NIR_vsb_TOA(i)
endif ! endif
write (6,*) 'SW down flux TOA = ',flx_SW_dwn_TOA(i),' W m-2'
write (6,*) 'SW up flux TOA = ',flx_SW_up_TOA(i),' W m-2'
write (6,*) 'SW net flux TOA = ',flx_SW_dwn_TOA(i)-flx_SW_up_TOA(i),' W m-2'
write (6,*) 'Vis down flux TOA = ',flx_vsb_dwn_TOA(i),' W m-2'
write (6,*) 'Vis up flux TOA = ',flx_vsb_up_TOA(i),' W m-2'
write (6,*) 'Vis net flux TOA = ',flx_vsb_dwn_TOA(i)-flx_vsb_up_TOA(i),' W m-2'
write (6,*) 'NIR down flux TOA = ',flx_NIR_dwn_TOA(i),' W m-2'
write (6,*) 'NIR up flux TOA = ',flx_NIR_up_TOA(i),' W m-2'
c write (6,*) 'NIR net flux TOA = ',fsnirtsq(i),' W m-2'
write (6,*) 'NIR net flux TOA = ',flx_NIR_dwn_TOA(i)-flx_NIR_up_TOA(i),' W m-2'
write (6,*) 'NIR net flux TOA N7 = ',fsnirt(i),' W m-2'
write (6,*) 'NIR net flux TOA N7 (clr) = ',fsnrtc(i),' W m-2'
write (6,*) ''
c Output SRB in terms of the "true" (unscaled) direct and diffuse fluxes
write (6,'(a)') 'Solar surface radiation budget components:'
if (flx_SW_dwn_sfc(i).le.0.0) then
write (6,*) 'SW alb sfc = Ill-defined'
else
write (6,*) 'SW alb sfc = ',alb_SW_sfc(i)
endif ! endif
if (flx_vsb_dwn_sfc(i).le.0.0) then
write (6,*) 'Vis alb sfc = Ill-defined'
else
write (6,*) 'Vis alb sfc = ',alb_vsb_sfc(i)
endif ! endif
if (flx_NIR_dwn_sfc(i).le.0.0) then
write (6,*) 'NIR alb sfc = Ill-defined'
else
write (6,*) 'NIR alb sfc = ',alb_NIR_sfc(i)
endif ! endif
write (6,*) 'SW down flux sfc = ',flx_SW_dwn_sfc(i),' W m-2'
write (6,*) 'SW down flux drc sfc = ',flx_SW_dwn_drc_sfc(i),' W m-2'
write (6,*) 'SW down flux dff sfc = ',flx_SW_dwn_dff_sfc(i),' W m-2'
if (flx_SW_dwn_drc_sfc(i).le.0.0) then
write (6,*) 'SW down flux dff/drc = Ill-defined'
else
write (6,*) 'SW down flux dff/drc = ',dff_drc_SW_sfc(i)
endif ! endif
write (6,*) 'Vis down flux sfc = ',flx_vsb_dwn_sfc(i),' W m-2'
write (6,*) 'Vis down flux drc sfc = ',flx_vsb_dwn_drc_sfc(i),' W m-2'
write (6,*) 'Vis down flux dff sfc = ',flx_vsb_dwn_dff_sfc(i),' W m-2'
if (flx_vsb_dwn_drc_sfc(i).le.0.0) then
write (6,*) 'Vis down flux dff/drc = Ill-defined'
else
write (6,*) 'Vis down flux dff/drc = ',dff_drc_vsb_sfc(i)
endif ! endif
write (6,*) 'NIR down flux sfc = ',flx_NIR_dwn_sfc(i),' W m-2'
write (6,*) 'NIR down flux drc sfc = ',flx_NIR_dwn_drc_sfc(i),' W m-2'
write (6,*) 'NIR down flux dff sfc = ',flx_NIR_dwn_dff_sfc(i),' W m-2'
if (flx_NIR_dwn_drc_sfc(i).le.0.0) then
write (6,*) 'NIR down flux dff/drc = Ill-defined'
else
write (6,*) 'NIR down flux dff/drc = ',dff_drc_NIR_sfc(i)
endif ! endif
write (6,*) ''
write (6,*) ''
#endif /* not CRM_SRB */
c The following output the SRB in terms of the scaled beam
c These may be useful to some investigators so we leave them here
c Do not use these numbers unless you are sure you know what they mean,
c i.e., unless you know how delta scaling repartitions downwelling flux.
c write (6,*) 'SW down flux sfc = ',sols(i)+soll(i)+solsd(i)+solld(i),' W m-2'
c write (6,*) 'SW down flux drc sfc = ',sols(i)+soll(i),' W m-2'
c write (6,*) 'SW down flux dff sfc = ',solsd(i)+solld(i),' W m-2'
c if (sols(i)+soll(i).le.0.0) then
c write (6,*) 'SW down flux dff/drc = Ill-defined'
c else
c write (6,*) 'SW down flux dff/drc = ',(solsd(i)+solld(i))/(sols(i)+soll(i))
c endif
c write (6,*) 'Vis down flux sfc = ',sols(i)+solsd(i),' W m-2'
c write (6,*) 'Vis down flux drc sfc = ',sols(i),' W m-2'
c write (6,*) 'Vis down flux dff sfc = ',solsd(i),' W m-2'
c if (sols(i).le.0.0) then
c write (6,*) 'Vis down flux dff/drc = Ill-defined'
c else
c write (6,*) 'Vis down flux dff/drc = ',solsd(i)/sols(i)
c endif
c write (6,*) 'NIR down flux sfc = ',soll(i)+solld(i),' W m-2'
c write (6,*) 'NIR down flux drc sfc = ',soll(i),' W m-2'
c write (6,*) 'NIR down flux dff sfc = ',solld(i),' W m-2'
c if (soll(i).le.0.0) then
c write (6,*) 'NIR down flux dff/drc = Ill-defined'
c else
c write (6,*) 'NIR down flux dff/drc = ',solld(i)/soll(i)
c endif
write (6,'(a)') 'Cloud microphysics:'
write (6,*) 'Level Pressure r_e lqd r_e ice Ice frc'
write (6,*) ' # mb um um frc'
do k=1,plev
write (6,'(i4,4x,f8.3,1x,3(f12.7,1x))')
$ k,pmidm1(i,k)/100.0,rel(i,k),rei(i,k),fice(i,k)
enddo ! end loop over levels
write (6,*) ''
#ifdef CRM_SRB
write (6,'(a)') 'SW Spectral Fluxes:'
write (6,*) 'Band Wvl Min Wvl Max Wvl Ctr TOA Dwn TOA Up Srf Dwn Srf Up'
write (6,*) ' # um um um W m-2 um-1 W m-2 um-1 W m-2 um-1 W m-2 um-1'
do k=1,bnd_nbr_SW
write (6,'(i4,1x,3(f8.4,1x),4(f12.7,1x))')
$ k,wvl_min(k)*1.0e6,wvl_max(k)*1.0e6,wvl_ctr(k)*1.0e6,
$ flx_bnd_dwn_TOA(i,k)*1.0e-6/wvl_dlt(k),flx_bnd_up_TOA(i,k)*1.0e-6/wvl_dlt(k),
$ flx_bnd_dwn_sfc(i,k)*1.0e-6/wvl_dlt(k),flx_bnd_up_sfc(i,k)*1.0e-6/wvl_dlt(k)
enddo ! end loop over bands
write (6,*) ''
write (6,'(a)') 'SW Scattering:'
write (6,*) 'Interface Pressure SW down SW direct SW diffuse SW dff/drc'
write (6,*) ' # mb W m-2 W m-2 W m-2 frc'
do k=1,plevp ! NB: plevp not plev
write (6,'(i4,8x,f8.3,1x,4(f12.7,1x))')
$ k,pintm1(i,k)/100.0,flx_SW_dwn(i,k),flx_SW_dwn_drc(i,k),
$ flx_SW_dwn_dff(i,k),dff_drc_SW(i,k)
enddo ! end loop over levels
write (6,*) ''
write (6,'(a)') 'SW Fluxes:'
write (6,*) 'Interface Pressure SW down SW up SW Net'
write (6,*) ' # mb W m-2 W m-2 W m-2'
do k=1,plevp ! NB: plevp not plev
write (6,'(i4,8x,f8.3,1x,3(f12.7,1x))')
$ k,pintm1(i,k)/100.0,flx_SW_dwn(i,k),flx_SW_up(i,k),
$ flx_SW_dwn(i,k)-flx_SW_up(i,k)
enddo ! end loop over levels
write (6,*) ''
write (6,'(a)') 'LW Fluxes:'
write (6,*) 'Interface Pressure LW down LW up LW Net'
write (6,*) ' # mb W m-2 W m-2 W m-2'
do k=1,plevp ! NB: plevp not plev
write (6,'(i4,8x,f8.3,1x,3(f12.7,1x))')
$ k,pintm1(i,k)/100.0,flx_LW_dwn(i,k),flx_LW_up(i,k),
$ flx_LW_up(i,k)-flx_LW_dwn(i,k)
enddo ! end loop over levels
write (6,*) ''
write (6,'(a)') 'Total SW+LW Fluxes:'
write (6,*) 'Interface Pressure Down Up Net'
write (6,*) ' # mb W m-2 W m-2 W m-2'
do k=1,plevp ! NB: plevp not plev
write (6,'(i4,8x,f8.3,1x,3(f12.7,1x))')
$ k,pintm1(i,k)/100.0,flx_SW_dwn(i,k)+flx_LW_dwn(i,k),
$ flx_SW_up(i,k)+flx_LW_up(i,k),
$ flx_SW_dwn(i,k)+flx_LW_dwn(i,k)-(flx_SW_up(i,k)+flx_LW_up(i,k))
enddo ! end loop over levels
write (6,*) ''
#endif /* not CRM_SRB */
write (6,'(a)') 'Heating rates:'
write (6,*) 'Level Pressure SW LW Net'
write (6,*) ' # mb K day-1 K day-1 K day-1'
do k=1,plev
write (6,'(i4,4x,f8.3,1x,3(f12.7,1x))')
$ k,pmidm1(i,k)/100.0,qrs(i,k)*86400.0,qrl(i,k)*86400.0,
$ 86400.0*(qrs(i,k)+qrl(i,k))
enddo ! end loop over levels
enddo ! end loop over longitude
c Write results to a netCDF file
call netcdf(
$ clat,
$ cld(1,1),
$ clon(1),
$ clwp(1,1),
$ coszrs(1),
$ effcld(1,1),
$ fice(1,1),
c
$ flns(1),
$ flnsc(1),
$ flnt(1),
$ flntc(1),
$ flwds(1),
c
$ fsnirt(1), ! [W m-2] Near-IR flux absorbed at TOA
$ fsnirtsq(1), ! [W m-2] Near-IR flux absorbed at TOA >= 0.7 microns
$ fsnrtc(1), ! [W m-2] Clear sky near-IR flux absorbed at TOA
c
$ fsns(1),
$ fsnsc(1),
$ fsnt(1),
$ fsntc(1),
$ loctim(1),
c
$ lwup(1),
$ pintm1(1,1),
$ pmidm1(1,1),
$ ps(1),
$ o3mmr(1,1),
$ o3vmr(1,1),
c
$ oro(1),
$ qm1(1,1),
$ qrl(1,1),
$ qrs(1,1),
$ rei(1,1),
c
$ rel(1,1),
$ snowh(1),
$ solin(1),
$ srfrad(1),
$ tg(1),
c
$ tm1(1,1),
$ ts(1))
write (6,'(a)') 'End CCM3 CRM'
stop
end
subroutine getdat( 1,1
$ aldif,
$ aldir,
$ asdif,
$ asdir,
$ clat,
$ cld,
$ clon,
$ clwp,
$ coslat,
$ loctim,
$ o3mmr,
$ o3vmr,
$ oro,
$ pilnm1,
$ pintm1,
$ pmidm1,
$ pmlnm1,
$ ps,
$ qm1,
$ snowh,
$ tg,
$ tm1,
$ ts)
c Purpose: Initialize column thermodynamic profile from external data
c Variables in getdat() have the same names they have in tphysbc(), where possible
c O3 mass mixing ratios are read in, but the model also requires O3 path lengths; they are computed here
c Cloud longwave emissivity is computed from cloud input (fraction and liquid water path), this is done here
#include <implicit.h>
c Parameters
#include <prgrid.h>
c Commons
#include <comtim.h> /* calday */
#include <crdcon.h> /* pie */
#include <comvmr.h> /* co2vmr, n2ovmr, ch4vmr, cfc11, cfc12 */
#include <comsol.h> /* eccen, obliq, mvelp, iyear_AD */
c Output arguments
real aldif(plond) ! Albedo: longwave, diffuse
real aldir(plond) ! Albedo: longwave, direct
real asdif(plond) ! Albedo: shortwave, diffuse
real asdir(plond) ! Albedo: shortwave, direct
real clat ! model latitude in radians
real cld(plond,plevp) ! cloud fraction
real clon(plon) ! Centered longitude (radians)
real clwp(plond,plev) ! cloud liquid water path (g/m**2)
real coslat ! cosine latitude
real loctim(plond) ! local time of solar computation
real o3mmr(plond,plev) ! o3 mass mixing ratio
real o3vmr(plond,plev) ! o3 volume mixing ratio
real oro(plond) ! Land surface flag
real pilnm1(plond,plevp) ! ln(pintm1)
real pintm1(plond,plevp) ! pressure at model interfaces
real pmidm1(plond,plev) ! pressure at model mid-levels
real pmlnm1(plond,plev) ! ln(pmidm1)
real ps(plond) ! model surface pressure field
real qm1(plond,plev) ! moisture field
real snowh(plond) ! snow depth (liquid water equivalent)
real tg(plond) ! surface (skin) temperature
real tm1(plond,plev) ! atmospheric temperature
real ts(plond) ! surface (air) temperature
c Local workspace
character*80 lbl ! Temporary space for labels
integer dbg_lvl ! Debugging level
integer i ! longitude index
integer k ! level index
integer lev(plev) ! [mb] Level index input
logical log_print ! Flag for status information in orb_params()
real lat_dgr ! [dgr] Latitude input
real lon_dgr ! [dgr] Longitude input
real rghnss(plond) ! surface roughness (obsolete)
real frctst(plond) ! fraction of surface with strong zenith dependent albedo (obsolete)
c CCM: physics/radinp()
real amd ! effective molecular weight of dry air (g/mol)
real amo ! molecular weight of ozone (g/mol)
data amd / 28.9644 /
data amo / 48.0000 /
real vmmr ! Factor for ozone volume mixing ratio
c Main Code
c Initialize some variables
dbg_lvl=0
i=1 ! Longitude index
c Read input data
read (5,'(a80)') lbl
if (dbg_lvl.gt.0) write (6,*) lbl
read (5,'(a80)') lbl
if (dbg_lvl.gt.0) write (6,*) lbl
read (5,'(a80)') lbl
if (dbg_lvl.gt.0) write (6,*) lbl
read (5,*) calday
c CCM: control/comtim.h: calday set in CCM: ccmlsm_share/calendr()
if (dbg_lvl.gt.0) write (6,*) calday,' = Julian day of year (1.5 = Noon, Jan 1; from 1 to 365)'
read (5,*) lat_dgr
if (dbg_lvl.gt.0) write (6,*) lat_dgr,' = Latitude (degrees North, from -90.0 to +90.0)'
#ifdef NEW_FMT
read (5,*) lon_dgr
if (dbg_lvl.gt.0) write (6,*) lon_dgr,' = Longitude (degrees East, from 0.0 to 360.0)'
#endif /* not NEW_FMT */
read (5,'(a80)') lbl
if (dbg_lvl.gt.0) write (6,*) lbl
do k=1,plev
read (5,*) lev(k),pmidm1(i,k),tm1(i,k),qm1(i,k),o3mmr(i,k),cld(i,k),clwp(i,k)
if (dbg_lvl.gt.0) write (6,'(1x,i3,1x,6(1pe10.3,1x))') k,pmidm1(i,k),tm1(i,k),qm1(i,k),o3mmr(i,k),cld(i,k),clwp(i,k)
c CCM: physics/cldfrc()
cld(i,k)=min(cld(i,k),0.999)
enddo ! end loop over lev
read (5,*) ps(i)
if (dbg_lvl.gt.0) write (6,*) ps(i),' = Surface pressure [mb]'
c Convert pressures from mb to pascals and define interface pressures:
ps(i)=ps(i)*100.0 ! mb --> Pa
do k=1,plev
pmidm1(i,k)=pmidm1(i,k)*100.0 ! mb --> Pa
pmlnm1(i,k)=log(pmidm1(i,k))
enddo ! end loop over lev
do k=1,plevp
if(k.eq.1) then
pintm1(i,k)=pmidm1(i,k)/2.0 ! Pa
else if (k.gt.1.and.k.le.plev) then
pintm1(i,k)=0.5*(pmidm1(i,k-1)+pmidm1(i,k)) ! Pa
else if (k.eq.plevp) then
pintm1(i,k)=ps(i) ! Pa
endif ! end if
pilnm1(i,k)=log(pintm1(i,k))
enddo ! end loop over levp
read (5,*) ts(i)
if (dbg_lvl.gt.0) write (6,*) ts(i),' = Surface air temperature [K]'
read (5,*) tg(i)
if (dbg_lvl.gt.0) write (6,*) tg(i),' = Surface skin temperature [K]'
read (5,*) oro(i)
if (dbg_lvl.gt.0) write (6,*) oro(i),' = Surface type (ORO) flag'
read (5,*) rghnss(i)
if (dbg_lvl.gt.0) write (6,*) rghnss(i),' = Surface aerodynamic roughness [m] (obsolete)'
read (5,*) snowh(i)
if (dbg_lvl.gt.0) write (6,*) snowh(i),' = Snow depth (liq. equiv.) [m]'
read (5,*) asdir(i)
if (dbg_lvl.gt.0) write (6,*) asdir(i),' = Albedo (Vis, direct)'
read (5,*) asdif(i)
if (dbg_lvl.gt.0) write (6,*) asdif(i),' = Albedo (Vis, diffuse)'
read (5,*) aldir(i)
if (dbg_lvl.gt.0) write (6,*) aldir(i),' = Albedo (NIR, direct)'
read (5,*) aldif(i)
if (dbg_lvl.gt.0) write (6,*) aldif(i),' = Albedo (NIR, diffuse)'
read (5,*) frctst(i)
if (dbg_lvl.gt.0) write (6,*) frctst(i),' = Fraction strong zenith angle dep. sfc. (obsolete)'
read (5,*) co2vmr ! CCM: physics/comvmr.h: co2vmr set in control/preset()
if (dbg_lvl.gt.0) write (6,*) co2vmr,' = CO2 volume mixing ratio [# #-1]'
read (5,*) n2ovmr ! CCM: physics/comvmr.h: n2ovmr set in control/preset()
if (dbg_lvl.gt.0) write (6,*) n2ovmr,' = N2O volume mixing ratio [# #-1]'
read (5,*) ch4vmr ! CCM: physics/comvmr.h: ch4vmr set in control/preset()
if (dbg_lvl.gt.0) write (6,*) ch4vmr,' = CH4 volume mixing ratio [# #-1]'
read (5,*) f11vmr ! CCM: physics/comvmr.h: co2vmr set in control/preset()
if (dbg_lvl.gt.0) write (6,*) f11vmr,' = CFC11 volume mixing ratio [# #-1]'
read (5,*) f12vmr ! CCM: physics/comvmr.h: f12vmr set in control/preset()
if (dbg_lvl.gt.0) write (6,*) f12vmr,' = CFC12 volume mixing ratio [# #-1]'
read (5,*) tauvis ! CCM: physics/comsol.h: tauvis set in control/preset()
if (dbg_lvl.gt.0) write (6,*) tauvis,' = Aerosol visible opt. depth'
read (5,*) scon ! CCM: physics/comsol.h: scon set in control/preset()
if (dbg_lvl.gt.0) write (6,*) scon,' = Solar constant [W m-2]'
read (5,*) iyear_AD ! CCM: physics/comsol.h: iyear_AD set in control/preset()
if (dbg_lvl.gt.0) write (6,*) iyear_AD,' = Year AD'
read (5,*) lon_dgr
if (dbg_lvl.gt.0) write (6,*) lon_dgr,' = Longitude (degrees East, from 0.0 to 360.0)'
c CCM: physics/comsol.h: scon set in control/preset()
scon=scon*1000.0 ! SI W m-2 = J m-2 s-1 --> CGS dyne cm-2 s-1
loctim(i)=(calday-aint(calday))*24.0
c CCM: control/crdcon.h: pie set in physics/radini()
pie=4.0*atan(1.0)
c CCM: control/commap.h: clat set in dynamics/eul/initcom()
clat=lat_dgr*pie/180.0
c CCM: control/commap.h: clon set in dynamics/eul/initcom()
clon(i)=lon_dgr*pie/180.0
c CCM: dynamics/eul/linemsbc()
coslat=cos(clat)
c Given iyear_AD, orb_params() computes physics/comsol.h: obliqr, lambm0 and mvelpp
c CCM: control/initext()
log_print=.false.
call orb_params( iyear_AD , eccen , obliq, mvelp,
$ obliqr , lambm0 , mvelpp, log_print )
c Convert ozone mass mixing ratio to ozone volume mixing ratio.
c CCM: physics/radinp.F (inverted version)
vmmr=amo/amd
do k=1,plev
do i=1,plon
c o3mmr(i,k)=vmmr*o3vmr(i,k)
o3vmr(i,k)=o3mmr(i,k)/vmmr
end do
end do
if (dbg_lvl.gt.0) write (6,*) 'End of profile data input'
return
end ! end getdat()
subroutine writeric(nabem,absems,lngbuf,nrow) 1
c Dummy routine for write of abs/ems data
c writeric() is called by radclw()
return
end
subroutine readric(nabem,absems,lngbuf,nrow) 1
c Dummy routine for read of abs/ems data
c readric() is called by radclw()
return
end ! end readric()
subroutine outfld(name,tmp ,plond,lat,hbuf) 18
c Dummy routine for history tape write
c outfld() is called by radctl()
return
end ! end outfld()
subroutine radozn(lat ,pmid ,o3vmr ) 1
c Dummy routine for ozone read
c radozn() is called from radctl()
c CRM gets o3vmr from getdat() instead
return
end ! end radozn()