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()