/wrfv2_fire/phys/module_ra_cam.F
FORTRAN Legacy | 7901 lines | 4745 code | 545 blank | 2611 comment | 157 complexity | 1d2649a1c0515f0f3dc8a5f07df2f95b MD5 | raw file
Possible License(s): AGPL-1.0
Large files files are truncated, but you can click here to view the full file
- MODULE module_ra_cam
- use module_ra_cam_support
- use module_cam_support, only: endrun
- implicit none
- !
- ! A. Slingo's data for cloud particle radiative properties (from 'A GCM
- ! Parameterization for the Shortwave Properties of Water Clouds' JAS
- ! vol. 46 may 1989 pp 1419-1427)
- !
- real(r8) abarl(4) ! A coefficient for extinction optical depth
- real(r8) bbarl(4) ! B coefficient for extinction optical depth
- real(r8) cbarl(4) ! C coefficient for single scat albedo
- real(r8) dbarl(4) ! D coefficient for single scat albedo
- real(r8) ebarl(4) ! E coefficient for asymmetry parameter
- real(r8) fbarl(4) ! F coefficient for asymmetry parameter
- save abarl, bbarl, cbarl, dbarl, ebarl, fbarl
- data abarl/ 2.817e-02, 2.682e-02,2.264e-02,1.281e-02/
- data bbarl/ 1.305 , 1.346 ,1.454 ,1.641 /
- data cbarl/-5.62e-08 ,-6.94e-06 ,4.64e-04 ,0.201 /
- data dbarl/ 1.63e-07 , 2.35e-05 ,1.24e-03 ,7.56e-03 /
- data ebarl/ 0.829 , 0.794 ,0.754 ,0.826 /
- data fbarl/ 2.482e-03, 4.226e-03,6.560e-03,4.353e-03/
- #if 0
- ! moved and changed to local variables into radcswmx for thread-safety, JM 20100217
- real(r8) abarli ! A coefficient for current spectral band
- real(r8) bbarli ! B coefficient for current spectral band
- real(r8) cbarli ! C coefficient for current spectral band
- real(r8) dbarli ! D coefficient for current spectral band
- real(r8) ebarli ! E coefficient for current spectral band
- real(r8) fbarli ! F coefficient for current spectral band
- #endif
- !
- ! Caution... A. Slingo recommends no less than 4.0 micro-meters nor
- ! greater than 20 micro-meters
- !
- ! ice water coefficients (Ebert and Curry,1992, JGR, 97, 3831-3836)
- !
- real(r8) abari(4) ! a coefficient for extinction optical depth
- real(r8) bbari(4) ! b coefficient for extinction optical depth
- real(r8) cbari(4) ! c coefficient for single scat albedo
- real(r8) dbari(4) ! d coefficient for single scat albedo
- real(r8) ebari(4) ! e coefficient for asymmetry parameter
- real(r8) fbari(4) ! f coefficient for asymmetry parameter
- save abari, bbari, cbari, dbari, ebari, fbari
- data abari/ 3.448e-03, 3.448e-03,3.448e-03,3.448e-03/
- data bbari/ 2.431 , 2.431 ,2.431 ,2.431 /
- data cbari/ 1.00e-05 , 1.10e-04 ,1.861e-02,.46658 /
- data dbari/ 0.0 , 1.405e-05,8.328e-04,2.05e-05 /
- data ebari/ 0.7661 , 0.7730 ,0.794 ,0.9595 /
- data fbari/ 5.851e-04, 5.665e-04,7.267e-04,1.076e-04/
- #if 0
- ! moved and changed to local variables into radcswmx for thread-safety, JM 20100217
- real(r8) abarii ! A coefficient for current spectral band
- real(r8) bbarii ! B coefficient for current spectral band
- real(r8) cbarii ! C coefficient for current spectral band
- real(r8) dbarii ! D coefficient for current spectral band
- real(r8) ebarii ! E coefficient for current spectral band
- real(r8) fbarii ! F coefficient for current spectral band
- #endif
- !
- real(r8) delta ! Pressure (in atm) for stratos. h2o limit
- real(r8) o2mmr ! O2 mass mixing ratio:
- save delta, o2mmr
- !
- ! UPDATE TO H2O NEAR-IR: Delta optimized for Hitran 2K and CKD 2.4
- !
- data delta / 0.0014257179260883 /
- !
- ! END UPDATE
- !
- data o2mmr / .23143 /
- ! Next series depends on spectral interval
- !
- real(r8) frcsol(nspint) ! Fraction of solar flux in spectral interval
- real(r8) wavmin(nspint) ! Min wavelength (micro-meters) of interval
- real(r8) wavmax(nspint) ! Max wavelength (micro-meters) of interval
- real(r8) raytau(nspint) ! Rayleigh scattering optical depth
- real(r8) abh2o(nspint) ! Absorption coefficiant for h2o (cm2/g)
- real(r8) abo3 (nspint) ! Absorption coefficiant for o3 (cm2/g)
- real(r8) abco2(nspint) ! Absorption coefficiant for co2 (cm2/g)
- real(r8) abo2 (nspint) ! Absorption coefficiant for o2 (cm2/g)
- real(r8) ph2o(nspint) ! Weight of h2o in spectral interval
- real(r8) pco2(nspint) ! Weight of co2 in spectral interval
- real(r8) po2 (nspint) ! Weight of o2 in spectral interval
- real(r8) nirwgt(nspint) ! Spectral Weights to simulate Nimbus-7 filter
- save frcsol ,wavmin ,wavmax ,raytau ,abh2o ,abo3 , &
- abco2 ,abo2 ,ph2o ,pco2 ,po2 ,nirwgt
- data frcsol / .001488, .001389, .001290, .001686, .002877, &
- .003869, .026336, .360739, .065392, .526861, &
- .526861, .526861, .526861, .526861, .526861, &
- .526861, .006239, .001834, .001834/
- !
- ! weight for 0.64 - 0.7 microns appropriate to clear skies over oceans
- !
- data nirwgt / 0.0, 0.0, 0.0, 0.0, 0.0, &
- 0.0, 0.0, 0.0, 0.320518, 1.0, 1.0, &
- 1.0, 1.0, 1.0, 1.0, 1.0, &
- 1.0, 1.0, 1.0 /
- data wavmin / .200, .245, .265, .275, .285, &
- .295, .305, .350, .640, .700, .701, &
- .701, .701, .701, .702, .702, &
- 2.630, 4.160, 4.160/
- data wavmax / .245, .265, .275, .285, .295, &
- .305, .350, .640, .700, 5.000, 5.000, &
- 5.000, 5.000, 5.000, 5.000, 5.000, &
- 2.860, 4.550, 4.550/
- !
- ! UPDATE TO H2O NEAR-IR: Rayleigh scattering optimized for Hitran 2K & CKD 2.4
- !
- real(r8) v_raytau_35
- real(r8) v_raytau_64
- real(r8) v_abo3_35
- real(r8) v_abo3_64
- parameter( &
- v_raytau_35 = 0.155208, &
- v_raytau_64 = 0.0392, &
- v_abo3_35 = 2.4058030e+01, &
- v_abo3_64 = 2.210e+01 &
- )
- data raytau / 4.020, 2.180, 1.700, 1.450, 1.250, &
- 1.085, 0.730, v_raytau_35, v_raytau_64, &
- 0.02899756, 0.01356763, 0.00537341, &
- 0.00228515, 0.00105028, 0.00046631, &
- 0.00025734, &
- .0001, .0001, .0001/
- !
- ! END UPDATE
- !
- !
- ! Absorption coefficients
- !
- !
- ! UPDATE TO H2O NEAR-IR: abh2o optimized for Hitran 2K and CKD 2.4
- !
- data abh2o / .000, .000, .000, .000, .000, &
- .000, .000, .000, .000, &
- 0.00256608, 0.06310504, 0.42287445, 2.45397941, &
- 11.20070807, 47.66091389, 240.19010243, &
- .000, .000, .000/
- !
- ! END UPDATE
- !
- data abo3 /5.370e+04, 13.080e+04, 9.292e+04, 4.530e+04, 1.616e+04, &
- 4.441e+03, 1.775e+02, v_abo3_35, v_abo3_64, .000, &
- .000, .000 , .000 , .000 , .000, &
- .000, .000 , .000 , .000 /
- data abco2 / .000, .000, .000, .000, .000, &
- .000, .000, .000, .000, .000, &
- .000, .000, .000, .000, .000, &
- .000, .094, .196, 1.963/
- data abo2 / .000, .000, .000, .000, .000, &
- .000, .000, .000,1.11e-05,6.69e-05, &
- .000, .000, .000, .000, .000, &
- .000, .000, .000, .000/
- !
- ! Spectral interval weights
- !
- data ph2o / .000, .000, .000, .000, .000, &
- .000, .000, .000, .000, .505, &
- .210, .120, .070, .048, .029, &
- .018, .000, .000, .000/
- data pco2 / .000, .000, .000, .000, .000, &
- .000, .000, .000, .000, .000, &
- .000, .000, .000, .000, .000, &
- .000, 1.000, .640, .360/
- data po2 / .000, .000, .000, .000, .000, &
- .000, .000, .000, 1.000, 1.000, &
- .000, .000, .000, .000, .000, &
- .000, .000, .000, .000/
- real(r8) amo ! Molecular weight of ozone (g/mol)
- save amo
- data amo / 48.0000 /
- contains
- subroutine camrad(RTHRATENLW,RTHRATENSW, &
- dolw,dosw, &
- SWUPT,SWUPTC,SWDNT,SWDNTC, &
- LWUPT,LWUPTC,LWDNT,LWDNTC, &
- SWUPB,SWUPBC,SWDNB,SWDNBC, &
- LWUPB,LWUPBC,LWDNB,LWDNBC, &
- swcf,lwcf,olr,cemiss,taucldc,taucldi,coszr, &
- GSW,GLW,XLAT,XLONG, &
- ALBEDO,t_phy,TSK,EMISS, &
- QV3D,QC3D,QR3D,QI3D,QS3D,QG3D, &
- ALSWVISDIR,ALSWVISDIF, & !ssib
- ALSWNIRDIR,ALSWNIRDIF, & !ssib
- SWVISDIR,SWVISDIF, & !ssib
- SWNIRDIR,SWNIRDIF, & !ssib
- sf_surface_physics, & !ssib
- F_QV,F_QC,F_QR,F_QI,F_QS,F_QG, &
- f_ice_phy,f_rain_phy, &
- p_phy,p8w,z,pi_phy,rho_phy,dz8w, &
- CLDFRA,XLAND,XICE,SNOW, &
- ozmixm,pin0,levsiz,num_months, &
- m_psp,m_psn,aerosolcp,aerosolcn,m_hybi0, &
- cam_abs_dim1, cam_abs_dim2, &
- paerlev,naer_c, &
- GMT,JULDAY,JULIAN,YR,DT,XTIME,DECLIN,SOLCON, &
- RADT,DEGRAD,n_cldadv, &
- abstot_3d, absnxt_3d, emstot_3d, &
- doabsems, &
- ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte )
- USE module_wrf_error
- USE module_state_description, ONLY : SSIBSCHEME !ssib
- !------------------------------------------------------------------
- IMPLICIT NONE
- !------------------------------------------------------------------
- INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- LOGICAL, INTENT(IN ) :: F_QV,F_QC,F_QR,F_QI,F_QS,F_QG
- LOGICAL, INTENT(INout) :: doabsems
- LOGICAL, INTENT(IN ) :: dolw,dosw
- INTEGER, INTENT(IN ) :: n_cldadv
- INTEGER, INTENT(IN ) :: JULDAY
- REAL, INTENT(IN ) :: JULIAN
- INTEGER, INTENT(IN ) :: YR
- REAL, INTENT(IN ) :: DT
- INTEGER, INTENT(IN ) :: levsiz, num_months
- INTEGER, INTENT(IN ) :: paerlev, naer_c
- INTEGER, INTENT(IN ) :: cam_abs_dim1, cam_abs_dim2
- REAL, INTENT(IN ) :: RADT,DEGRAD, &
- XTIME,DECLIN,SOLCON,GMT
- !
- !
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
- INTENT(IN ) :: P_PHY, &
- P8W, &
- Z, &
- pi_PHY, &
- rho_PHY, &
- dz8w, &
- T_PHY, &
- QV3D, &
- QC3D, &
- QR3D, &
- QI3D, &
- QS3D, &
- QG3D, &
- CLDFRA
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
- INTENT(INOUT) :: RTHRATENLW, &
- RTHRATENSW
- !
- REAL, DIMENSION( ims:ime, jms:jme ), &
- INTENT(IN ) :: XLAT, &
- XLONG, &
- XLAND, &
- XICE, &
- SNOW, &
- EMISS, &
- TSK, &
- ALBEDO
- REAL, DIMENSION( ims:ime, levsiz, jms:jme, num_months ), &
- INTENT(IN ) :: OZMIXM
- REAL, DIMENSION(levsiz), INTENT(IN ) :: PIN0
- REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN ) :: m_psp,m_psn
- REAL, DIMENSION(paerlev), intent(in) :: m_hybi0
- REAL, DIMENSION( ims:ime, paerlev, jms:jme, naer_c ), &
- INTENT(IN ) :: aerosolcp, aerosolcn
- !
- REAL, DIMENSION( ims:ime, jms:jme ), &
- INTENT(INOUT) :: GSW, GLW
- !---------SSiB variables (fds 06/2010)----------------
- REAL, DIMENSION( ims:ime, jms:jme ), &
- INTENT(IN) :: ALSWVISDIR, &
- ALSWVISDIF, &
- ALSWNIRDIR, &
- ALSWNIRDIF
- REAL, DIMENSION( ims:ime, jms:jme ), &
- INTENT(OUT) :: SWVISDIR, &
- SWVISDIF, &
- SWNIRDIR, &
- SWNIRDIF
- INTEGER, INTENT(IN) :: sf_surface_physics
- !--------------------------------------
- ! saving arrays for doabsems reduction of radiation calcs
- REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim2 , jms:jme ), &
- INTENT(INOUT) :: abstot_3d
- REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim1 , jms:jme ), &
- INTENT(INOUT) :: absnxt_3d
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
- INTENT(INOUT) :: emstot_3d
- ! Added outputs of total and clearsky fluxes etc
- ! Note that k=1 refers to the half level below the model lowest level (Sfc)
- ! k=kme refers to the half level above the model highest level (TOA)
- !
- ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
- ! INTENT(INOUT) :: swup, &
- ! swupclear, &
- ! swdn, &
- ! swdnclear, &
- ! lwup, &
- ! lwupclear, &
- ! lwdn, &
- ! lwdnclear
- REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
- SWUPT,SWUPTC,SWDNT,SWDNTC, &
- LWUPT,LWUPTC,LWDNT,LWDNTC, &
- SWUPB,SWUPBC,SWDNB,SWDNBC, &
- LWUPB,LWUPBC,LWDNB,LWDNBC
- REAL, DIMENSION( ims:ime, jms:jme ), &
- INTENT(INOUT) :: swcf, &
- lwcf, &
- olr, &
- coszr
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
- INTENT(OUT ) :: cemiss, & ! cloud emissivity for isccp
- taucldc, & ! cloud water optical depth for isccp
- taucldi ! cloud ice optical depth for isccp
- !
- !
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
- INTENT(IN ) :: &
- F_ICE_PHY, &
- F_RAIN_PHY
- ! LOCAL VARIABLES
-
- INTEGER :: lchnk, ncol, pcols, pver, pverp, pverr, pverrp
- INTEGER :: pcnst, pnats, ppcnst, i, j, k, ii, kk, kk1, m, n
- integer :: begchunk, endchunk
- integer :: nyrm, nyrp
- real(r8) doymodel, doydatam, doydatap, deltat, fact1, fact2
- REAL :: XT24, TLOCTM, HRANG, XXLAT, oldXT24
-
- real(r8), DIMENSION( 1:ite-its+1 ) :: coszrs, landfrac, landm, snowh, icefrac, lwups
- real(r8), DIMENSION( 1:ite-its+1 ) :: asdir, asdif, aldir, aldif, ps
- real(r8), DIMENSION( 1:ite-its+1, 1:kte-kts+1 ) :: cld, pmid, lnpmid, pdel, zm, t
- real(r8), DIMENSION( 1:ite-its+1, 1:kte-kts+2 ) :: pint, lnpint
- real(r8), DIMENSION( 1:ite-its+1, 1:kte-kts+1, n_cldadv) :: q
- ! real(r8), DIMENSION( 1:kte-kts+1 ) :: hypm ! reference pressures at midpoints
- ! real(r8), DIMENSION( 1:kte-kts+2 ) :: hypi ! reference pressures at interfaces
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: cicewp ! in-cloud cloud ice water path
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: cliqwp ! in-cloud cloud liquid water path
- real(r8), dimension( 1:ite-its+1, 0:kte-kts+1 ) :: tauxcl ! cloud water optical depth
- real(r8), dimension( 1:ite-its+1, 0:kte-kts+1 ) :: tauxci ! cloud ice optical depth
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: emis ! cloud emissivity
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: rel ! effective drop radius (microns)
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: rei ! ice effective drop size (microns)
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: pmxrgn ! Maximum values of pressure for each
- integer , dimension( 1:ite-its+1 ) :: nmxrgn ! Number of maximally overlapped regions
- real(r8), dimension( 1:ite-its+1 ) :: fsns ! Surface absorbed solar flux
- real(r8), dimension( 1:ite-its+1 ) :: fsnt ! Net column abs solar flux at model top
- real(r8), dimension( 1:ite-its+1 ) :: flns ! Srf longwave cooling (up-down) flux
- real(r8), dimension( 1:ite-its+1 ) :: flnt ! Net outgoing lw flux at model top
- ! Added outputs of total and clearsky fluxes etc
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fsup ! Upward total sky solar
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fsupc ! Upward clear sky solar
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fsdn ! Downward total sky solar
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fsdnc ! Downward clear sky solar
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: flup ! Upward total sky longwave
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: flupc ! Upward clear sky longwave
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fldn ! Downward total sky longwave
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fldnc ! Downward clear sky longwave
- real(r8), dimension( 1:ite-its+1 ) :: swcftoa ! Top of the atmosphere solar cloud forcing
- real(r8), dimension( 1:ite-its+1 ) :: lwcftoa ! Top of the atmosphere longwave cloud forcing
- real(r8), dimension( 1:ite-its+1 ) :: olrtoa ! Top of the atmosphere outgoing longwave
- !
- real(r8), dimension( 1:ite-its+1 ) :: sols ! Downward solar rad onto surface (sw direct)
- real(r8), dimension( 1:ite-its+1 ) :: soll ! Downward solar rad onto surface (lw direct)
- real(r8), dimension( 1:ite-its+1 ) :: solsd ! Downward solar rad onto surface (sw diffuse)
- real(r8), dimension( 1:ite-its+1 ) :: solld ! Downward solar rad onto surface (lw diffuse)
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: qrs ! Solar heating rate
- real(r8), dimension( 1:ite-its+1 ) :: fsds ! Flux Shortwave Downwelling Surface
- real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: qrl ! Longwave cooling rate
- real(r8), dimension( 1:ite-its+1 ) :: flwds ! Surface down longwave flux
- real(r8), dimension( 1:ite-its+1, levsiz, num_months ) :: ozmixmj ! monthly ozone mixing ratio
- real(r8), dimension( 1:ite-its+1, levsiz ) :: ozmix ! ozone mixing ratio (time interpolated)
- real(r8), dimension(levsiz) :: pin ! ozone pressure level
- real(r8), dimension(1:ite-its+1) :: m_psjp,m_psjn ! MATCH surface pressure
- real(r8), dimension( 1:ite-its+1, paerlev, naer_c ) :: aerosoljp ! monthly aerosol concentrations
- real(r8), dimension( 1:ite-its+1, paerlev, naer_c ) :: aerosoljn ! monthly aerosol concentrations
- real(r8), dimension(paerlev) :: m_hybi
- real(r8), dimension(1:ite-its+1 ) :: clat ! latitude in radians for columns
- real(r8), dimension(its:ite,kts:kte+1,kts:kte+1) :: abstot ! Total absorptivity
- real(r8), dimension(its:ite,kts:kte,4) :: absnxt ! Total nearest layer absorptivity
- real(r8), dimension(its:ite,kts:kte+1) :: emstot ! Total emissivity
- CHARACTER(LEN=256) :: msgstr
- #if !defined(MAC_KLUDGE)
- lchnk = 1
- begchunk = ims
- endchunk = ime
- ncol = ite - its + 1
- pcols= ite - its + 1
- pver = kte - kts + 1
- pverp= pver + 1
- pverr = kte - kts + 1
- pverrp= pverr + 1
- ! number of advected constituents and non-advected constituents (including water vapor)
- ppcnst = n_cldadv
- ! number of non-advected constituents
- pnats = 0
- pcnst = ppcnst-pnats
- ! check the # species defined for the input climatology and naer
- ! if(naer_c.ne.naer) then
- ! WRITE( wrf_err_message , * ) 'naer_c ne naer ', naer_c, naer
- if(naer_c.ne.naer_all) then
- WRITE( wrf_err_message , * ) 'naer_c-1 ne naer_all ', naer_c, naer_all
- CALL wrf_error_fatal ( wrf_err_message )
- endif
- ! update CO2 volume mixing ratio (co2vmr)
-
- ! determine time interpolation factors, check sanity
- ! of interpolation factors to within 32-bit roundoff
- ! assume that day of year is 1 for all input data
- !
- nyrm = yr - yrdata(1) + 1
- nyrp = nyrm + 1
- doymodel = yr*365. + julian
- doydatam = yrdata(nyrm)*365. + 1.
- doydatap = yrdata(nyrp)*365. + 1.
- deltat = doydatap - doydatam
- fact1 = (doydatap - doymodel)/deltat
- fact2 = (doymodel - doydatam)/deltat
- co2vmr = (co2(nyrm)*fact1 + co2(nyrp)*fact2)*1.e-06
- co2mmr=co2vmr*mwco2/mwdry
- !
- !===================================================
- ! Radiation computations
- !===================================================
- do k=1,levsiz
- pin(k)=pin0(k)
- enddo
- do k=1,paerlev
- m_hybi(k)=m_hybi0(k)
- enddo
- ! check for uninitialized arrays
- if(abstot_3d(its,kts,kts,jts) .eq. 0.0 .and. .not.doabsems .and. dolw)then
- CALL wrf_debug(0, 'camrad lw: CAUTION: re-calculating abstot, absnxt, emstot on restart')
- doabsems = .true.
- endif
- do j =jts,jte
- !
- ! Cosine solar zenith angle for current time step
- !
- ! call zenith (calday, clat, clon, coszrs, ncol)
- do i = its,ite
- ii = i - its + 1
- ! XT24 is the fractional part of simulation days plus half of RADT expressed in
- ! units of minutes
- ! JULIAN is in days
- ! RADT is in minutes
- XT24=MOD(XTIME+RADT*0.5,1440.)
- TLOCTM=GMT+XT24/60.+XLONG(I,J)/15.
- HRANG=15.*(TLOCTM-12.)*DEGRAD
- XXLAT=XLAT(I,J)*DEGRAD
- clat(ii)=xxlat
- coszrs(II)=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG)
- enddo
- ! moist variables
- do k = kts,kte
- kk = kte - k + kts
- do i = its,ite
- ii = i - its + 1
- ! convert to specific humidity
- q(ii,kk,1) = max(1.e-10,qv3d(i,k,j)/(1.+qv3d(i,k,j)))
- IF ( F_QI .and. F_QC .and. F_QS ) THEN
- q(ii,kk,ixcldliq) = max(0.,qc3d(i,k,j)/(1.+qv3d(i,k,j)))
- q(ii,kk,ixcldice) = max(0.,(qi3d(i,k,j)+qs3d(i,k,j))/(1.+qv3d(i,k,j)))
- ELSE IF ( F_QC .and. F_QR ) THEN
- ! Warm rain or simple ice
- q(ii,kk,ixcldliq) = 0.
- q(ii,kk,ixcldice) = 0.
- if(t_phy(i,k,j).gt.273.15)q(ii,kk,ixcldliq) = max(0.,qc3d(i,k,j)/(1.+qv3d(i,k,j)))
- if(t_phy(i,k,j).le.273.15)q(ii,kk,ixcldice) = max(0.,qc3d(i,k,j)/(1.+qv3d(i,k,j)))
- ELSE IF ( F_QC .and. F_QS ) THEN
- ! For Ferrier (note that currently Ferrier has QI, so this section will not be used)
- q(ii,kk,ixcldice) = max(0.,qc3d(i,k,j)/(1.+qv3d(i,k,j))*f_ice_phy(i,k,j))
- q(ii,kk,ixcldliq) = max(0.,qc3d(i,k,j)/(1.+qv3d(i,k,j))*(1.-f_ice_phy(i,k,j))*(1.-f_rain_phy(i,k,j)))
- ELSE
- q(ii,kk,ixcldliq) = 0.
- q(ii,kk,ixcldice) = 0.
- ENDIF
- cld(ii,kk) = CLDFRA(I,K,J)
- enddo
- enddo
- do i = its,ite
- ii = i - its + 1
- landfrac(ii) = 2.-XLAND(I,J)
- landm(ii) = landfrac(ii)
- snowh(ii) = 0.001*SNOW(I,J)
- icefrac(ii) = XICE(I,J)
- enddo
- do m=1,num_months-1
- do k=1,levsiz
- do i = its,ite
- ii = i - its + 1
- ozmixmj(ii,k,m) = ozmixm(i,k,j,m+1)
- enddo
- enddo
- enddo
- do i = its,ite
- ii = i - its + 1
- m_psjp(ii) = m_psp(i,j)
- m_psjn(ii) = m_psn(i,j)
- enddo
- do n=1,naer_c
- do k=1,paerlev
- do i = its,ite
- ii = i - its + 1
- aerosoljp(ii,k,n) = aerosolcp(i,k,j,n)
- aerosoljn(ii,k,n) = aerosolcn(i,k,j,n)
- enddo
- enddo
- enddo
- !
- ! Complete radiation calculations
- !
- do i = its,ite
- ii = i - its + 1
- lwups(ii) = stebol*EMISS(I,J)*TSK(I,J)**4
- enddo
- do k = kts,kte+1
- kk = kte - k + kts + 1
- do i = its,ite
- ii = i - its + 1
- pint(ii,kk) = p8w(i,k,j)
- if(k.eq.kts)ps(ii)=pint(ii,kk)
- lnpint(ii,kk) = log(pint(ii,kk))
- enddo
- enddo
- if(.not.doabsems .and. dolw)then
- ! do kk = kts,kte+1
- do kk = 1,cam_abs_dim2
- do kk1 = kts,kte+1
- do i = its,ite
- abstot(i,kk1,kk) = abstot_3d(i,kk1,kk,j)
- enddo
- enddo
- enddo
- ! do kk = 1,4
- do kk = 1,cam_abs_dim1
- do kk1 = kts,kte
- do i = its,ite
- absnxt(i,kk1,kk) = absnxt_3d(i,kk1,kk,j)
- enddo
- enddo
- enddo
- do kk = kts,kte+1
- do i = its,ite
- emstot(i,kk) = emstot_3d(i,kk,j)
- enddo
- enddo
- endif
- do k = kts,kte
- kk = kte - k + kts
- do i = its,ite
- ii = i - its + 1
- pmid(ii,kk) = p_phy(i,k,j)
- lnpmid(ii,kk) = log(pmid(ii,kk))
- lnpint(ii,kk) = log(pint(ii,kk))
- pdel(ii,kk) = pint(ii,kk+1) - pint(ii,kk)
- t(ii,kk) = t_phy(i,k,j)
- zm(ii,kk) = z(i,k,j)
- enddo
- enddo
- ! Compute cloud water/ice paths and optical properties for input to radiation
- call param_cldoptics_calc(ncol, pcols, pver, pverp, pverr, pverrp, ppcnst, q, cld, landfrac, landm,icefrac, &
- pdel, t, ps, pmid, pint, cicewp, cliqwp, emis, rel, rei, pmxrgn, nmxrgn, snowh)
- !-----fds (06/2010)----------------------------
- SELECT CASE(sf_surface_physics)
- CASE (SSIBSCHEME)
- if (xtime .gt. 1.0) then
- ! call wrf_message("using SSiB albedoes for land points")
- do i = its,ite
- ii = i - its + 1
- if (xland(i,j).lt.1.5) then !land points only
- asdir(ii) = ALSWVISDIR(i,j) ! SSiB visdir albedo
- asdif(ii) = ALSWVISDIF(i,j) ! SSiB visdif albedo
- aldir(ii) = ALSWNIRDIR(i,j) ! SSiB nirdir albedo
- aldif(ii) = ALSWNIRDIF(i,j) ! SSiB nirdif albedo
- else
- asdir(ii) = albedo(i,j)
- asdif(ii) = albedo(i,j)
- aldir(ii) = albedo(i,j)
- aldif(ii) = albedo(i,j)
- endif
- enddo
- else
- do i = its,ite
- ii = i - its + 1
- asdir(ii) = albedo(i,j)
- asdif(ii) = albedo(i,j)
- aldir(ii) = albedo(i,j)
- aldif(ii) = albedo(i,j)
- enddo
- endif
- CASE DEFAULT
- do i = its,ite
- ii = i - its + 1
- ! use same albedo for direct and diffuse
- ! change this when separate values are provided
- asdir(ii) = albedo(i,j)
- asdif(ii) = albedo(i,j)
- aldir(ii) = albedo(i,j)
- aldif(ii) = albedo(i,j)
- enddo
- END SELECT
- !-----------------------------------------------
- ! WRF allocate space here (not needed if oznini is called)
- ! allocate (ozmix(pcols,levsiz,begchunk:endchunk)) ! This line from oznini.F90
- call radctl (j,lchnk, ncol, pcols, pver, pverp, pverr, pverrp, ppcnst, pcnst, lwups, emis, pmid, &
- pint, lnpmid, lnpint, pdel, t, q, &
- cld, cicewp, cliqwp, tauxcl, tauxci, coszrs, clat, asdir, asdif, &
- aldir, aldif, solcon, GMT,JULDAY,JULIAN,DT,XTIME, &
- pin, ozmixmj, ozmix, levsiz, num_months, &
- m_psjp,m_psjn, aerosoljp, aerosoljn, m_hybi, paerlev, naer_c, pmxrgn, nmxrgn, &
- dolw, dosw, doabsems, abstot, absnxt, emstot, &
- fsup, fsupc, fsdn, fsdnc, flup, flupc, fldn, fldnc, swcftoa, lwcftoa, olrtoa, &
- fsns, fsnt ,flns ,flnt , &
- qrs, qrl, flwds, rel, rei, &
- sols, soll, solsd, solld, &
- landfrac, zm, fsds)
- do k = kts,kte
- kk = kte - k + kts
- do i = its,ite
- ii = i - its + 1
- if(dolw)RTHRATENLW(I,K,J) = 1.e4*qrl(ii,kk)/(cpair*pi_phy(i,k,j))
- if(dosw)RTHRATENSW(I,K,J) = 1.e4*qrs(ii,kk)/(cpair*pi_phy(i,k,j))
- cemiss(i,k,j) = emis(ii,kk)
- taucldc(i,k,j) = tauxcl(ii,kk)
- taucldi(i,k,j) = tauxci(ii,kk)
- enddo
- enddo
- if(doabsems .and. dolw)then
- ! do kk = kts,kte+1
- do kk = 1,cam_abs_dim2
- do kk1 = kts,kte+1
- do i = its,ite
- abstot_3d(i,kk1,kk,j) = abstot(i,kk1,kk)
- enddo
- enddo
- enddo
- ! do kk = 1,4
- do kk = 1,cam_abs_dim1
- do kk1 = kts,kte
- do i = its,ite
- absnxt_3d(i,kk1,kk,j) = absnxt(i,kk1,kk)
- enddo
- enddo
- enddo
- do kk = kts,kte+1
- do i = its,ite
- emstot_3d(i,kk,j) = emstot(i,kk)
- enddo
- enddo
- endif
- IF(PRESENT(SWUPT))THEN
- if(dosw)then
- ! Added shortwave and longwave upward/downward total and clear sky fluxes
- do k = kts,kte+1
- kk = kte +1 - k + kts
- do i = its,ite
- ii = i - its + 1
- ! swup(i,k,j) = fsup(ii,kk)
- ! swupclear(i,k,j) = fsupc(ii,kk)
- ! swdn(i,k,j) = fsdn(ii,kk)
- ! swdnclear(i,k,j) = fsdnc(ii,kk)
- if(k.eq.kte+1)then
- swupt(i,j) = fsup(ii,kk)
- swuptc(i,j) = fsupc(ii,kk)
- swdnt(i,j) = fsdn(ii,kk)
- swdntc(i,j) = fsdnc(ii,kk)
- endif
- if(k.eq.kts)then
- swupb(i,j) = fsup(ii,kk)
- swupbc(i,j) = fsupc(ii,kk)
- swdnb(i,j) = fsdn(ii,kk)
- swdnbc(i,j) = fsdnc(ii,kk)
- endif
- ! if(i.eq.30.and.j.eq.30) then
- ! print 1234, 'short ', i,ii,k,kk,fsup(ii,kk),fsupc(ii,kk),fsdn(ii,kk),fsdnc(ii,kk)
- ! 1234 format (a6,4i4,4f10.3)
- ! endif
- enddo
- enddo
- endif
- if(dolw)then
- ! Added shortwave and longwave upward/downward total and clear sky fluxes
- do k = kts,kte+1
- kk = kte +1 - k + kts
- do i = its,ite
- ii = i - its + 1
- ! lwup(i,k,j) = flup(ii,kk)
- ! lwupclear(i,k,j) = flupc(ii,kk)
- ! lwdn(i,k,j) = fldn(ii,kk)
- ! lwdnclear(i,k,j) = fldnc(ii,kk)
- if(k.eq.kte+1)then
- lwupt(i,j) = flup(ii,kk)
- lwuptc(i,j) = flupc(ii,kk)
- lwdnt(i,j) = fldn(ii,kk)
- lwdntc(i,j) = fldnc(ii,kk)
- endif
- if(k.eq.kts)then
- lwupb(i,j) = flup(ii,kk)
- lwupbc(i,j) = flupc(ii,kk)
- lwdnb(i,j) = fldn(ii,kk)
- lwdnbc(i,j) = fldnc(ii,kk)
- endif
- ! if(i.eq.30.and.j.eq.30) then
- ! print 1234, 'long ', i,ii,k,kk,flup(ii,kk),flupc(ii,kk),fldn(ii,kk),fldnc(ii,kk)
- ! 1234 format (a6,4i4,4f10.3)
- ! endif
- enddo
- enddo
- endif
- ENDIF
- do i = its,ite
- ii = i - its + 1
- ! Added shortwave and longwave cloud forcing at TOA and surface
- if(dolw)then
- GLW(I,J) = flwds(ii)
- lwcf(i,j) = lwcftoa(ii)
- olr(i,j) = olrtoa(ii)
- endif
- if(dosw)then
- GSW(I,J) = fsns(ii)
- swcf(i,j) = swcftoa(ii)
- coszr(i,j) = coszrs(ii)
- endif
- enddo
- !-------fds (06/2010)---------
- SELECT CASE(sf_surface_physics)
- CASE (SSIBSCHEME)
- ! call wrf_message("CAM using ssib albedo2")
- if(dosw)then
- do i = its,ite
- ii = i - its + 1
- SWVISDIR(I,J) = sols(ii) !SSiB
- SWVISDIF(I,J) = solsd(ii) !SSiB
- SWNIRDIR(I,J) = soll(ii) !SSiB
- SWNIRDIF(I,J) = solld(ii) !SSiB
- enddo
- endif
- END SELECT
- !-----------------------------
- enddo ! j-loop
- #endif
- end subroutine camrad
- !====================================================================
- SUBROUTINE camradinit( &
- R_D,R_V,CP,G,STBOLT,EP_2,shalf,pptop, &
- ozmixm,pin,levsiz,XLAT,num_months, &
- m_psp,m_psn,m_hybi,aerosolcp,aerosolcn, &
- paerlev,naer_c, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte )
- USE module_wrf_error
- USE module_state_description
- !USE module_configure
- !--------------------------------------------------------------------
- IMPLICIT NONE
- !--------------------------------------------------------------------
- INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte
- REAL, intent(in) :: pptop
- REAL, INTENT(IN) :: R_D,R_V,CP,G,STBOLT,EP_2
- REAL, DIMENSION( kms:kme ) :: shalf
- INTEGER, INTENT(IN ) :: levsiz, num_months
- INTEGER, INTENT(IN ) :: paerlev, naer_c
- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: XLAT
- REAL, DIMENSION( ims:ime, levsiz, jms:jme, num_months ), &
- INTENT(INOUT ) :: OZMIXM
- REAL, DIMENSION(levsiz), INTENT(INOUT ) :: PIN
- REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT ) :: m_psp,m_psn
- REAL, DIMENSION(paerlev), INTENT(INOUT ) :: m_hybi
- REAL, DIMENSION( ims:ime, paerlev, jms:jme, naer_c ), &
- INTENT(INOUT) :: aerosolcp,aerosolcn
- REAL(r8) :: pstd
- REAL(r8) :: rh2o, cpair
- ! These were made allocatable 20090612 to save static memory allocation. JM
- IF ( .NOT. ALLOCATED( ksul ) ) ALLOCATE( ksul( nrh, nspint ) )
- IF ( .NOT. ALLOCATED( wsul ) ) ALLOCATE( wsul( nrh, nspint ) )
- IF ( .NOT. ALLOCATED( gsul ) ) ALLOCATE( gsul( nrh, nspint ) )
- IF ( .NOT. ALLOCATED( ksslt ) ) ALLOCATE( ksslt( nrh, nspint ) )
- IF ( .NOT. ALLOCATED( wsslt ) ) ALLOCATE( wsslt( nrh, nspint ) )
- IF ( .NOT. ALLOCATED( gsslt ) ) ALLOCATE( gsslt( nrh, nspint ) )
- IF ( .NOT. ALLOCATED( kcphil ) ) ALLOCATE( kcphil( nrh, nspint ) )
- IF ( .NOT. ALLOCATED( wcphil ) ) ALLOCATE( wcphil( nrh, nspint ) )
- IF ( .NOT. ALLOCATED( gcphil ) ) ALLOCATE( gcphil( nrh, nspint ) )
- IF ( .NOT. ALLOCATED(ah2onw ) ) ALLOCATE( ah2onw(n_p, n_tp, n_u, n_te, n_rh) )
- IF ( .NOT. ALLOCATED(eh2onw ) ) ALLOCATE( eh2onw(n_p, n_tp, n_u, n_te, n_rh) )
- IF ( .NOT. ALLOCATED(ah2ow ) ) ALLOCATE( ah2ow(n_p, n_tp, n_u, n_te, n_rh) )
- IF ( .NOT. ALLOCATED(cn_ah2ow) ) ALLOCATE( cn_ah2ow(n_p, n_tp, n_u, n_te, n_rh) )
- IF ( .NOT. ALLOCATED(cn_eh2ow) ) ALLOCATE( cn_eh2ow(n_p, n_tp, n_u, n_te, n_rh) )
- IF ( .NOT. ALLOCATED(ln_ah2ow) ) ALLOCATE( ln_ah2ow(n_p, n_tp, n_u, n_te, n_rh) )
- IF ( .NOT. ALLOCATED(ln_eh2ow) ) ALLOCATE( ln_eh2ow(n_p, n_tp, n_u, n_te, n_rh) )
- #if !defined(MAC_KLUDGE)
- ozncyc = .true.
- indirect = .true.
- ixcldliq = 2
- ixcldice = 3
- #if (NMM_CORE != 1)
- ! aerosol array is not in the NMM Registry
- ! since CAM radiation not available to NMM (yet)
- ! so this is blocked out to enable CAM compilation with NMM
- idxSUL = P_SUL
- idxSSLT = P_SSLT
- idxDUSTfirst = P_DUST1
- idxOCPHO = P_OCPHO
- idxCARBONfirst = P_OCPHO
- idxBCPHO = P_BCPHO
- idxOCPHI = P_OCPHI
- idxBCPHI = P_BCPHI
- idxBG = P_BG
- idxVOLC = P_VOLC
- #endif
- pstd = 101325.0
- ! from physconst module
- mwdry = 28.966 ! molecular weight dry air ~ kg/kmole (shr_const_mwdair)
- mwco2 = 44. ! molecular weight co2
- mwh2o = 18.016 ! molecular weight water vapor (shr_const_mwwv)
- mwch4 = 16. ! molecular weight ch4
- mwn2o = 44. ! molecular weight n2o
- mwf11 = 136. ! molecular weight cfc11
- mwf12 = 120. ! molecular weight cfc12
- cappa = R_D/CP
- rair = R_D
- tmelt = 273.16 ! freezing T of fresh water ~ K
- r_universal = 6.02214e26 * STBOLT ! Universal gas constant ~ J/K/kmole
- latvap = 2.501e6 ! latent heat of evaporation ~ J/kg
- latice = 3.336e5 ! latent heat of fusion ~ J/kg
- zvir = R_V/R_D - 1.
- rh2o = R_V
- cpair = CP
- !
- epsqs = EP_2
- CALL radini(G, CP, EP_2, STBOLT, pstd*10.0 )
- CALL esinti(epsqs ,latvap ,latice ,rh2o ,cpair ,tmelt )
- CALL oznini(ozmixm,pin,levsiz,num_months,XLAT, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte)
- CALL aerosol_init(m_psp,m_psn,m_hybi,aerosolcp,aerosolcn,paerlev,naer_c,shalf,pptop, &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte)
- #endif
- END SUBROUTINE camradinit
- #if !defined(MAC_KLUDGE)
- subroutine oznint(julday,julian,dt,gmt,xtime,ozmixmj,ozmix,levsiz,num_months,pcols)
- IMPLICIT NONE
- INTEGER, INTENT(IN ) :: levsiz, num_months,pcols
- REAL(r8), DIMENSION( pcols, levsiz, num_months ), &
- INTENT(IN ) :: ozmixmj
- REAL, INTENT(IN ) :: XTIME,GMT
- INTEGER, INTENT(IN ) :: JULDAY
- REAL, INTENT(IN ) :: JULIAN
- REAL, INTENT(IN ) :: DT
- REAL(r8), DIMENSION( pcols, levsiz ), &
- INTENT(OUT ) :: ozmix
- !Local
- REAL(r8) :: intJULIAN
- integer :: np1,np,nm,m,k,i
- integer :: IJUL
- integer, dimension(12) :: date_oz
- data date_oz/16, 45, 75, 105, 136, 166, 197, 228, 258, 289, 319, 350/
- real(r8) :: cdayozp, cdayozm
- real(r8) :: fact1, fact2
- logical :: finddate
- CHARACTER(LEN=256) :: msgstr
- ! JULIAN starts from 0.0 at 0Z on 1 Jan.
- intJULIAN = JULIAN + 1.0_r8 ! offset by one day
- ! jan 1st 00z is julian=1.0 here
- IJUL=INT(intJULIAN)
- ! Note that following will drift.
- ! Need to use actual month/day info to compute julian.
- intJULIAN=intJULIAN-FLOAT(IJUL)
- IJUL=MOD(IJUL,365)
- IF(IJUL.EQ.0)IJUL=365
- intJULIAN=intJULIAN+IJUL
- np1=1
- finddate=.false.
- ! do m=1,num_months
- do m=1,12
- if(date_oz(m).gt.intjulian.and..not.finddate) then
- np1=m
- finddate=.true.
- endif
- enddo
- cdayozp=date_oz(np1)
- if(np1.gt.1) then
- cdayozm=date_oz(np1-1)
- np=np1
- nm=np-1
- else
- cdayozm=date_oz(12)
- np=np1
- nm=12
- endif
- call getfactors(ozncyc,np1, cdayozm, cdayozp,intjulian, &
- fact1, fact2)
- !
- ! Time interpolation.
- !
- do k=1,levsiz
- do i=1,pcols
- ozmix(i,k) = ozmixmj(i,k,nm)*fact1 + ozmixmj(i,k,np)*fact2
- end do
- end do
- END subroutine oznint
- subroutine get_aerosol(c, julday, julian, dt, gmt, xtime, m_psp, m_psn, aerosoljp, &
- aerosoljn, m_hybi, paerlev, naer_c, pint, pcols, pver, pverp, pverr, pverrp, AEROSOLt, scale)
- !------------------------------------------------------------------
- !
- ! Input:
- ! time at which aerosol mmrs are needed (get_curr_calday())
- ! chunk index
- ! CAM's vertical grid (pint)
- !
- ! Output:
- ! values for Aerosol Mass Mixing Ratios at specified time
- ! on vertical grid specified by CAM (AEROSOLt)
- !
- ! Method:
- ! first determine which indexs of aerosols are the bounding data sets
- ! interpolate both onto vertical grid aerm(),aerp().
- ! from those two, interpolate in time.
- !
- !------------------------------------------------------------------
- ! use volcanicmass, only: get_volcanic_mass
- ! use timeinterp, only: getfactors
- !
- ! aerosol fields interpolated to current time step
- ! on pressure levels of this time step.
- ! these should be made read-only for other modules
- ! Is allocation done correctly here?
- !
- integer, intent(in) :: c ! Chunk Id.
- integer, intent(in) :: paerlev, naer_c, pcols, pver, pverp, pverr, pverrp
- real(r8), intent(in) :: pint(pcols,pverp) ! midpoint pres.
- real(r8), intent(in) :: scale(naer_all) ! scale each aerosol by this amount
- REAL, INTENT(IN ) :: XTIME,GMT
- INTEGER, INTENT(IN ) :: JULDAY
- REAL, INTENT(IN ) :: JULIAN
- REAL, INTENT(IN ) :: DT
- real(r8), intent(in ) :: m_psp(pcols),m_psn(pcols) ! Match surface pressure
- real(r8), intent(in ) :: aerosoljp(pcols,paerlev,naer_c)
- real(r8), intent(in ) :: aerosoljn(pcols,paerlev,naer_c)
- real(r8), intent(in ) :: m_hybi(paerlev)
- real(r8), intent(out) :: AEROSOLt(pcols, pver, naer_all) ! aerosols
- !
- ! Local workspace
- !
- real(r8) caldayloc ! calendar day of current timestep
- real(r8) fact1, fact2 ! time interpolation factors
- integer :: nm = 1 ! index to prv month in array. init to 1 and toggle between 1 and 2
- integer :: np = 2 ! index to nxt month in array. init to 2 and toggle between 1 and 2
- integer :: mo_nxt = bigint ! index to nxt month in file
- integer :: mo_prv ! index to previous month
- real(r8) :: cdaym = inf ! calendar day of prv month
- real(r8) :: cdayp = inf ! calendar day of next month
- real(r8) :: Mid(12) ! Days into year for mid month date
- data Mid/16.5, 46.0, 75.5, 106.0, 136.5, 167.0, 197.5, 228.5, 259.0, 289.5, 320.0, 350.5 /
- integer i, k, j ! spatial indices
- integer m ! constituent index
- integer lats(pcols),lons(pcols) ! latitude and longitudes of column
- integer ncol ! number of columns
- INTEGER IJUL
- REAL(r8) intJULIAN
- real(r8) speciesmin(naer) ! minimal value for each species
- !
- ! values before current time step "the minus month"
- ! aerosolm(pcols,pver) is value of preceeding month's aerosol mmr
- ! aerosolp(pcols,pver) is value of next month's aerosol mmr
- ! (think minus and plus or values to left and right of point to be interpolated)
- !
- real(r8) AEROSOLm(pcols,pver,naer) ! aerosol mmr from MATCH in column at previous (minus) month
- !
- ! values beyond (or at) current time step "the plus month"
- !
- real(r8) AEROSOLp(pcols,pver,naer) ! aerosol mmr from MATCH in column at next (plus) month
- CHARACTER(LEN=256) :: msgstr
- ! JULIAN starts from 0.0 at 0Z on 1 Jan.
- intJULIAN = JULIAN + 1.0_r8 ! offset by one day
- ! jan 1st 00z is julian=1.0 here
- IJUL=INT(intJULIAN)
- ! Note that following will drift.
- ! Need to use actual month/day info to compute julian.
- intJULIAN=intJULIAN-FLOAT(IJUL)
- IJUL=MOD(IJUL,365)
- IF(IJUL.EQ.0)IJUL=365
- caldayloc=intJULIAN+IJUL
- if (caldayloc < Mid(1)) then
- mo_prv = 12
- mo_nxt = 1
- else if (caldayloc >= Mid(12)) then
- mo_prv = 12
- mo_nxt = 1
- else
- do i = 2 , 12
- if (caldayloc < Mid(i)) then
- mo_prv = i-1
- mo_nxt = i
- exit
- end if
- end do
- end if
- !
- ! Set initial calendar day values
- !
- cdaym = Mid(mo_prv)
- cdayp = Mid(mo_nxt)
- !
- ! Determine time interpolation factors. 1st arg says we are cycling 1 year of data
- !
- call getfactors (.true., mo_nxt, cdaym, cdayp, caldayloc, &
- fact1, fact2)
- !
- ! interpolate (prv and nxt month) bounding datasets onto cam vertical grid.
- ! compute mass mixing ratios on CAMS's pressure coordinate
- ! for both the "minus" and "plus" months
- !
- ! ncol = get_ncols_p(c)
- ncol = pcols
- ! call vert_interpolate (M_ps_cam_col(1,c,nm), pint, nm, AEROSOLm, ncol, c)
- ! call vert_interpolate (M_ps_cam_col(1,c,np), pint, np, AEROSOLp, ncol, c)
- call vert_interpolate (m_psp, aerosoljp, m_hybi, paerlev, naer_c, pint, nm, AEROSOLm, pcols, pver, pverp, ncol, c)
- call vert_interpolate (m_psn, aerosoljn, m_hybi, paerlev, naer_c, pint, np, AEROSOLp, pcols, pver, pverp, ncol, c)
- !
- ! Time interpolate.
- !
- do m=1,naer
- do k=1,pver
- do i=1,ncol
- AEROSOLt(i,k,m) = AEROSOLm(i,k,m)*fact1 + AEROSOLp(i,k,m)*fact2
- end do
- end do
- end do
- ! do i=1,ncol
- ! Match_ps_chunk(i,c) = m_ps(i,nm)*fact1 + m_ps(i,np)*fact2
- ! end do
- !
- ! get background aerosol (tuning) field
- !
- call background (c, ncol, pint, pcols, pverr, pverrp, AEROSOLt(:, :, idxBG))
- !
- ! find volcanic aerosol masses
- !
- ! if (strat_volcanic) then
- ! call get_volcanic_mass(c, AEROSOLt(:,:,idxVOLC))
- ! else
- AEROSOLt(:,:,idxVOLC) = 0._r8
- ! endif
- !
- ! exit if mmr is negative (we have previously set
- ! cumulative mass to be a decreasing function.)
- !
- speciesmin(:) = 0. ! speciesmin(m) = 0 is minimum mmr for each species
- do m=1,naer
- do k=1,pver
- do i=1,ncol
- if (AEROSOLt(i, k, m) < speciesmin(m)) then
- write(6,*) 'AEROSOL_INTERPOLATE: negative mass mixing ratio, exiting'
- write(6,*) 'm, column, pver',m, i, k ,AEROSOLt(i, k, m)
- call endrun ()
- end if
- end do
- end do
- end do
- !
- ! scale any AEROSOLS as required
- !
- call scale_aerosols (AEROSOLt, pcols, pver, ncol, c, scale)
- return
- end subroutine get_aerosol
- subroutine aerosol_indirect(ncol,lchnk,pcols,pver,ppcnst,landfrac,pmid,t,qm1,cld,zm,rel)
- !--------------------------------------------------------------
- ! Compute effect of sulfate on effective liquid water radius
- ! Method of Martin et. al.
- !--------------------------------------------------------------
- ! use constituents, only: ppcnst, cnst_get_ind
- ! use history, only: outfld
- !#include <comctl.h>
- integer, intent(in) :: ncol ! number of atmospheric columns
- integer, intent(in) :: lchnk ! chunk identifier
- integer, intent(in) :: pcols,pver,ppcnst
- real(r8), intent(in) :: landfrac(pcols) ! land fraction
- real(r8), intent(in) :: pmid(pcols,pver) ! Model level pressures
- real(r8), intent(in) :: t(pcols,pver) ! Model level temperatures
- real(r8), intent(in) :: qm1(pcols,pver,ppcnst) ! Specific humidity and tracers
- real(r8), intent(in) :: cld(pcols,pver) ! Fractional cloud cover
- real(r8), intent(in) :: zm(pcols,pver) ! Height of midpoints (above surface)
- real(r8), intent(in) :: rel(pcols,pver) ! liquid effective drop size (microns)
- !
- ! local variables
- !
- real(r8) locrhoair(pcols,pver) ! dry air density [kg/m^3 ]
- real(r8) lwcwat(pcols,pver) ! in-cloud liquid water path [kg/m^3 ]
- real(r8) sulfmix(pcols,pver) ! sulfate mass mixing ratio [kg/kg ]
- real(r8) so4mass(pcols,pver) ! sulfate mass concentration [g/cm^3 ]
- real(r8) Aso4(pcols,pver) ! sulfate # concentration [#/cm^3 ]
- real(r8) Ntot(pcols,pver) ! ccn # concentration [#/cm^3 ]
- real(r8) relmod(pcols,pver) ! effective radius [microns]
- real(r8) wrel(pcols,pver) ! weighted effective radius [microns]
- real(r8) wlwc(pcols,pver) ! weighted liq. water content [kg/m^3 ]
- real(r8) cldfrq(pcols,pver) ! frequency of occurance of...
- ! ! clouds (cld => 0.01) [fraction]
- real(r8) locPi ! my piece of the pi
- real(r8) Rdryair ! gas constant of dry air [J/deg/kg]
- real(r8) rhowat ! density of water [kg/m^3 ]
- re…
Large files files are truncated, but you can click here to view the full file