PageRenderTime 71ms CodeModel.GetById 33ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_ra_cam.F

http://github.com/jbeezley/wrf-fire
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

  1. MODULE module_ra_cam
  2. use module_ra_cam_support
  3. use module_cam_support, only: endrun
  4. implicit none
  5. !
  6. ! A. Slingo's data for cloud particle radiative properties (from 'A GCM
  7. ! Parameterization for the Shortwave Properties of Water Clouds' JAS
  8. ! vol. 46 may 1989 pp 1419-1427)
  9. !
  10. real(r8) abarl(4) ! A coefficient for extinction optical depth
  11. real(r8) bbarl(4) ! B coefficient for extinction optical depth
  12. real(r8) cbarl(4) ! C coefficient for single scat albedo
  13. real(r8) dbarl(4) ! D coefficient for single scat albedo
  14. real(r8) ebarl(4) ! E coefficient for asymmetry parameter
  15. real(r8) fbarl(4) ! F coefficient for asymmetry parameter
  16. save abarl, bbarl, cbarl, dbarl, ebarl, fbarl
  17. data abarl/ 2.817e-02, 2.682e-02,2.264e-02,1.281e-02/
  18. data bbarl/ 1.305 , 1.346 ,1.454 ,1.641 /
  19. data cbarl/-5.62e-08 ,-6.94e-06 ,4.64e-04 ,0.201 /
  20. data dbarl/ 1.63e-07 , 2.35e-05 ,1.24e-03 ,7.56e-03 /
  21. data ebarl/ 0.829 , 0.794 ,0.754 ,0.826 /
  22. data fbarl/ 2.482e-03, 4.226e-03,6.560e-03,4.353e-03/
  23. #if 0
  24. ! moved and changed to local variables into radcswmx for thread-safety, JM 20100217
  25. real(r8) abarli ! A coefficient for current spectral band
  26. real(r8) bbarli ! B coefficient for current spectral band
  27. real(r8) cbarli ! C coefficient for current spectral band
  28. real(r8) dbarli ! D coefficient for current spectral band
  29. real(r8) ebarli ! E coefficient for current spectral band
  30. real(r8) fbarli ! F coefficient for current spectral band
  31. #endif
  32. !
  33. ! Caution... A. Slingo recommends no less than 4.0 micro-meters nor
  34. ! greater than 20 micro-meters
  35. !
  36. ! ice water coefficients (Ebert and Curry,1992, JGR, 97, 3831-3836)
  37. !
  38. real(r8) abari(4) ! a coefficient for extinction optical depth
  39. real(r8) bbari(4) ! b coefficient for extinction optical depth
  40. real(r8) cbari(4) ! c coefficient for single scat albedo
  41. real(r8) dbari(4) ! d coefficient for single scat albedo
  42. real(r8) ebari(4) ! e coefficient for asymmetry parameter
  43. real(r8) fbari(4) ! f coefficient for asymmetry parameter
  44. save abari, bbari, cbari, dbari, ebari, fbari
  45. data abari/ 3.448e-03, 3.448e-03,3.448e-03,3.448e-03/
  46. data bbari/ 2.431 , 2.431 ,2.431 ,2.431 /
  47. data cbari/ 1.00e-05 , 1.10e-04 ,1.861e-02,.46658 /
  48. data dbari/ 0.0 , 1.405e-05,8.328e-04,2.05e-05 /
  49. data ebari/ 0.7661 , 0.7730 ,0.794 ,0.9595 /
  50. data fbari/ 5.851e-04, 5.665e-04,7.267e-04,1.076e-04/
  51. #if 0
  52. ! moved and changed to local variables into radcswmx for thread-safety, JM 20100217
  53. real(r8) abarii ! A coefficient for current spectral band
  54. real(r8) bbarii ! B coefficient for current spectral band
  55. real(r8) cbarii ! C coefficient for current spectral band
  56. real(r8) dbarii ! D coefficient for current spectral band
  57. real(r8) ebarii ! E coefficient for current spectral band
  58. real(r8) fbarii ! F coefficient for current spectral band
  59. #endif
  60. !
  61. real(r8) delta ! Pressure (in atm) for stratos. h2o limit
  62. real(r8) o2mmr ! O2 mass mixing ratio:
  63. save delta, o2mmr
  64. !
  65. ! UPDATE TO H2O NEAR-IR: Delta optimized for Hitran 2K and CKD 2.4
  66. !
  67. data delta / 0.0014257179260883 /
  68. !
  69. ! END UPDATE
  70. !
  71. data o2mmr / .23143 /
  72. ! Next series depends on spectral interval
  73. !
  74. real(r8) frcsol(nspint) ! Fraction of solar flux in spectral interval
  75. real(r8) wavmin(nspint) ! Min wavelength (micro-meters) of interval
  76. real(r8) wavmax(nspint) ! Max wavelength (micro-meters) of interval
  77. real(r8) raytau(nspint) ! Rayleigh scattering optical depth
  78. real(r8) abh2o(nspint) ! Absorption coefficiant for h2o (cm2/g)
  79. real(r8) abo3 (nspint) ! Absorption coefficiant for o3 (cm2/g)
  80. real(r8) abco2(nspint) ! Absorption coefficiant for co2 (cm2/g)
  81. real(r8) abo2 (nspint) ! Absorption coefficiant for o2 (cm2/g)
  82. real(r8) ph2o(nspint) ! Weight of h2o in spectral interval
  83. real(r8) pco2(nspint) ! Weight of co2 in spectral interval
  84. real(r8) po2 (nspint) ! Weight of o2 in spectral interval
  85. real(r8) nirwgt(nspint) ! Spectral Weights to simulate Nimbus-7 filter
  86. save frcsol ,wavmin ,wavmax ,raytau ,abh2o ,abo3 , &
  87. abco2 ,abo2 ,ph2o ,pco2 ,po2 ,nirwgt
  88. data frcsol / .001488, .001389, .001290, .001686, .002877, &
  89. .003869, .026336, .360739, .065392, .526861, &
  90. .526861, .526861, .526861, .526861, .526861, &
  91. .526861, .006239, .001834, .001834/
  92. !
  93. ! weight for 0.64 - 0.7 microns appropriate to clear skies over oceans
  94. !
  95. data nirwgt / 0.0, 0.0, 0.0, 0.0, 0.0, &
  96. 0.0, 0.0, 0.0, 0.320518, 1.0, 1.0, &
  97. 1.0, 1.0, 1.0, 1.0, 1.0, &
  98. 1.0, 1.0, 1.0 /
  99. data wavmin / .200, .245, .265, .275, .285, &
  100. .295, .305, .350, .640, .700, .701, &
  101. .701, .701, .701, .702, .702, &
  102. 2.630, 4.160, 4.160/
  103. data wavmax / .245, .265, .275, .285, .295, &
  104. .305, .350, .640, .700, 5.000, 5.000, &
  105. 5.000, 5.000, 5.000, 5.000, 5.000, &
  106. 2.860, 4.550, 4.550/
  107. !
  108. ! UPDATE TO H2O NEAR-IR: Rayleigh scattering optimized for Hitran 2K & CKD 2.4
  109. !
  110. real(r8) v_raytau_35
  111. real(r8) v_raytau_64
  112. real(r8) v_abo3_35
  113. real(r8) v_abo3_64
  114. parameter( &
  115. v_raytau_35 = 0.155208, &
  116. v_raytau_64 = 0.0392, &
  117. v_abo3_35 = 2.4058030e+01, &
  118. v_abo3_64 = 2.210e+01 &
  119. )
  120. data raytau / 4.020, 2.180, 1.700, 1.450, 1.250, &
  121. 1.085, 0.730, v_raytau_35, v_raytau_64, &
  122. 0.02899756, 0.01356763, 0.00537341, &
  123. 0.00228515, 0.00105028, 0.00046631, &
  124. 0.00025734, &
  125. .0001, .0001, .0001/
  126. !
  127. ! END UPDATE
  128. !
  129. !
  130. ! Absorption coefficients
  131. !
  132. !
  133. ! UPDATE TO H2O NEAR-IR: abh2o optimized for Hitran 2K and CKD 2.4
  134. !
  135. data abh2o / .000, .000, .000, .000, .000, &
  136. .000, .000, .000, .000, &
  137. 0.00256608, 0.06310504, 0.42287445, 2.45397941, &
  138. 11.20070807, 47.66091389, 240.19010243, &
  139. .000, .000, .000/
  140. !
  141. ! END UPDATE
  142. !
  143. data abo3 /5.370e+04, 13.080e+04, 9.292e+04, 4.530e+04, 1.616e+04, &
  144. 4.441e+03, 1.775e+02, v_abo3_35, v_abo3_64, .000, &
  145. .000, .000 , .000 , .000 , .000, &
  146. .000, .000 , .000 , .000 /
  147. data abco2 / .000, .000, .000, .000, .000, &
  148. .000, .000, .000, .000, .000, &
  149. .000, .000, .000, .000, .000, &
  150. .000, .094, .196, 1.963/
  151. data abo2 / .000, .000, .000, .000, .000, &
  152. .000, .000, .000,1.11e-05,6.69e-05, &
  153. .000, .000, .000, .000, .000, &
  154. .000, .000, .000, .000/
  155. !
  156. ! Spectral interval weights
  157. !
  158. data ph2o / .000, .000, .000, .000, .000, &
  159. .000, .000, .000, .000, .505, &
  160. .210, .120, .070, .048, .029, &
  161. .018, .000, .000, .000/
  162. data pco2 / .000, .000, .000, .000, .000, &
  163. .000, .000, .000, .000, .000, &
  164. .000, .000, .000, .000, .000, &
  165. .000, 1.000, .640, .360/
  166. data po2 / .000, .000, .000, .000, .000, &
  167. .000, .000, .000, 1.000, 1.000, &
  168. .000, .000, .000, .000, .000, &
  169. .000, .000, .000, .000/
  170. real(r8) amo ! Molecular weight of ozone (g/mol)
  171. save amo
  172. data amo / 48.0000 /
  173. contains
  174. subroutine camrad(RTHRATENLW,RTHRATENSW, &
  175. dolw,dosw, &
  176. SWUPT,SWUPTC,SWDNT,SWDNTC, &
  177. LWUPT,LWUPTC,LWDNT,LWDNTC, &
  178. SWUPB,SWUPBC,SWDNB,SWDNBC, &
  179. LWUPB,LWUPBC,LWDNB,LWDNBC, &
  180. swcf,lwcf,olr,cemiss,taucldc,taucldi,coszr, &
  181. GSW,GLW,XLAT,XLONG, &
  182. ALBEDO,t_phy,TSK,EMISS, &
  183. QV3D,QC3D,QR3D,QI3D,QS3D,QG3D, &
  184. ALSWVISDIR,ALSWVISDIF, & !ssib
  185. ALSWNIRDIR,ALSWNIRDIF, & !ssib
  186. SWVISDIR,SWVISDIF, & !ssib
  187. SWNIRDIR,SWNIRDIF, & !ssib
  188. sf_surface_physics, & !ssib
  189. F_QV,F_QC,F_QR,F_QI,F_QS,F_QG, &
  190. f_ice_phy,f_rain_phy, &
  191. p_phy,p8w,z,pi_phy,rho_phy,dz8w, &
  192. CLDFRA,XLAND,XICE,SNOW, &
  193. ozmixm,pin0,levsiz,num_months, &
  194. m_psp,m_psn,aerosolcp,aerosolcn,m_hybi0, &
  195. cam_abs_dim1, cam_abs_dim2, &
  196. paerlev,naer_c, &
  197. GMT,JULDAY,JULIAN,YR,DT,XTIME,DECLIN,SOLCON, &
  198. RADT,DEGRAD,n_cldadv, &
  199. abstot_3d, absnxt_3d, emstot_3d, &
  200. doabsems, &
  201. ids,ide, jds,jde, kds,kde, &
  202. ims,ime, jms,jme, kms,kme, &
  203. its,ite, jts,jte, kts,kte )
  204. USE module_wrf_error
  205. USE module_state_description, ONLY : SSIBSCHEME !ssib
  206. !------------------------------------------------------------------
  207. IMPLICIT NONE
  208. !------------------------------------------------------------------
  209. INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
  210. ims,ime, jms,jme, kms,kme, &
  211. its,ite, jts,jte, kts,kte
  212. LOGICAL, INTENT(IN ) :: F_QV,F_QC,F_QR,F_QI,F_QS,F_QG
  213. LOGICAL, INTENT(INout) :: doabsems
  214. LOGICAL, INTENT(IN ) :: dolw,dosw
  215. INTEGER, INTENT(IN ) :: n_cldadv
  216. INTEGER, INTENT(IN ) :: JULDAY
  217. REAL, INTENT(IN ) :: JULIAN
  218. INTEGER, INTENT(IN ) :: YR
  219. REAL, INTENT(IN ) :: DT
  220. INTEGER, INTENT(IN ) :: levsiz, num_months
  221. INTEGER, INTENT(IN ) :: paerlev, naer_c
  222. INTEGER, INTENT(IN ) :: cam_abs_dim1, cam_abs_dim2
  223. REAL, INTENT(IN ) :: RADT,DEGRAD, &
  224. XTIME,DECLIN,SOLCON,GMT
  225. !
  226. !
  227. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
  228. INTENT(IN ) :: P_PHY, &
  229. P8W, &
  230. Z, &
  231. pi_PHY, &
  232. rho_PHY, &
  233. dz8w, &
  234. T_PHY, &
  235. QV3D, &
  236. QC3D, &
  237. QR3D, &
  238. QI3D, &
  239. QS3D, &
  240. QG3D, &
  241. CLDFRA
  242. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
  243. INTENT(INOUT) :: RTHRATENLW, &
  244. RTHRATENSW
  245. !
  246. REAL, DIMENSION( ims:ime, jms:jme ), &
  247. INTENT(IN ) :: XLAT, &
  248. XLONG, &
  249. XLAND, &
  250. XICE, &
  251. SNOW, &
  252. EMISS, &
  253. TSK, &
  254. ALBEDO
  255. REAL, DIMENSION( ims:ime, levsiz, jms:jme, num_months ), &
  256. INTENT(IN ) :: OZMIXM
  257. REAL, DIMENSION(levsiz), INTENT(IN ) :: PIN0
  258. REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN ) :: m_psp,m_psn
  259. REAL, DIMENSION(paerlev), intent(in) :: m_hybi0
  260. REAL, DIMENSION( ims:ime, paerlev, jms:jme, naer_c ), &
  261. INTENT(IN ) :: aerosolcp, aerosolcn
  262. !
  263. REAL, DIMENSION( ims:ime, jms:jme ), &
  264. INTENT(INOUT) :: GSW, GLW
  265. !---------SSiB variables (fds 06/2010)----------------
  266. REAL, DIMENSION( ims:ime, jms:jme ), &
  267. INTENT(IN) :: ALSWVISDIR, &
  268. ALSWVISDIF, &
  269. ALSWNIRDIR, &
  270. ALSWNIRDIF
  271. REAL, DIMENSION( ims:ime, jms:jme ), &
  272. INTENT(OUT) :: SWVISDIR, &
  273. SWVISDIF, &
  274. SWNIRDIR, &
  275. SWNIRDIF
  276. INTEGER, INTENT(IN) :: sf_surface_physics
  277. !--------------------------------------
  278. ! saving arrays for doabsems reduction of radiation calcs
  279. REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim2 , jms:jme ), &
  280. INTENT(INOUT) :: abstot_3d
  281. REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim1 , jms:jme ), &
  282. INTENT(INOUT) :: absnxt_3d
  283. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
  284. INTENT(INOUT) :: emstot_3d
  285. ! Added outputs of total and clearsky fluxes etc
  286. ! Note that k=1 refers to the half level below the model lowest level (Sfc)
  287. ! k=kme refers to the half level above the model highest level (TOA)
  288. !
  289. ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
  290. ! INTENT(INOUT) :: swup, &
  291. ! swupclear, &
  292. ! swdn, &
  293. ! swdnclear, &
  294. ! lwup, &
  295. ! lwupclear, &
  296. ! lwdn, &
  297. ! lwdnclear
  298. REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
  299. SWUPT,SWUPTC,SWDNT,SWDNTC, &
  300. LWUPT,LWUPTC,LWDNT,LWDNTC, &
  301. SWUPB,SWUPBC,SWDNB,SWDNBC, &
  302. LWUPB,LWUPBC,LWDNB,LWDNBC
  303. REAL, DIMENSION( ims:ime, jms:jme ), &
  304. INTENT(INOUT) :: swcf, &
  305. lwcf, &
  306. olr, &
  307. coszr
  308. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
  309. INTENT(OUT ) :: cemiss, & ! cloud emissivity for isccp
  310. taucldc, & ! cloud water optical depth for isccp
  311. taucldi ! cloud ice optical depth for isccp
  312. !
  313. !
  314. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
  315. INTENT(IN ) :: &
  316. F_ICE_PHY, &
  317. F_RAIN_PHY
  318. ! LOCAL VARIABLES
  319. INTEGER :: lchnk, ncol, pcols, pver, pverp, pverr, pverrp
  320. INTEGER :: pcnst, pnats, ppcnst, i, j, k, ii, kk, kk1, m, n
  321. integer :: begchunk, endchunk
  322. integer :: nyrm, nyrp
  323. real(r8) doymodel, doydatam, doydatap, deltat, fact1, fact2
  324. REAL :: XT24, TLOCTM, HRANG, XXLAT, oldXT24
  325. real(r8), DIMENSION( 1:ite-its+1 ) :: coszrs, landfrac, landm, snowh, icefrac, lwups
  326. real(r8), DIMENSION( 1:ite-its+1 ) :: asdir, asdif, aldir, aldif, ps
  327. real(r8), DIMENSION( 1:ite-its+1, 1:kte-kts+1 ) :: cld, pmid, lnpmid, pdel, zm, t
  328. real(r8), DIMENSION( 1:ite-its+1, 1:kte-kts+2 ) :: pint, lnpint
  329. real(r8), DIMENSION( 1:ite-its+1, 1:kte-kts+1, n_cldadv) :: q
  330. ! real(r8), DIMENSION( 1:kte-kts+1 ) :: hypm ! reference pressures at midpoints
  331. ! real(r8), DIMENSION( 1:kte-kts+2 ) :: hypi ! reference pressures at interfaces
  332. real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: cicewp ! in-cloud cloud ice water path
  333. real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: cliqwp ! in-cloud cloud liquid water path
  334. real(r8), dimension( 1:ite-its+1, 0:kte-kts+1 ) :: tauxcl ! cloud water optical depth
  335. real(r8), dimension( 1:ite-its+1, 0:kte-kts+1 ) :: tauxci ! cloud ice optical depth
  336. real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: emis ! cloud emissivity
  337. real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: rel ! effective drop radius (microns)
  338. real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: rei ! ice effective drop size (microns)
  339. real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: pmxrgn ! Maximum values of pressure for each
  340. integer , dimension( 1:ite-its+1 ) :: nmxrgn ! Number of maximally overlapped regions
  341. real(r8), dimension( 1:ite-its+1 ) :: fsns ! Surface absorbed solar flux
  342. real(r8), dimension( 1:ite-its+1 ) :: fsnt ! Net column abs solar flux at model top
  343. real(r8), dimension( 1:ite-its+1 ) :: flns ! Srf longwave cooling (up-down) flux
  344. real(r8), dimension( 1:ite-its+1 ) :: flnt ! Net outgoing lw flux at model top
  345. ! Added outputs of total and clearsky fluxes etc
  346. real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fsup ! Upward total sky solar
  347. real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fsupc ! Upward clear sky solar
  348. real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fsdn ! Downward total sky solar
  349. real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fsdnc ! Downward clear sky solar
  350. real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: flup ! Upward total sky longwave
  351. real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: flupc ! Upward clear sky longwave
  352. real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fldn ! Downward total sky longwave
  353. real(r8), dimension( 1:ite-its+1, 1:kte-kts+2 ) :: fldnc ! Downward clear sky longwave
  354. real(r8), dimension( 1:ite-its+1 ) :: swcftoa ! Top of the atmosphere solar cloud forcing
  355. real(r8), dimension( 1:ite-its+1 ) :: lwcftoa ! Top of the atmosphere longwave cloud forcing
  356. real(r8), dimension( 1:ite-its+1 ) :: olrtoa ! Top of the atmosphere outgoing longwave
  357. !
  358. real(r8), dimension( 1:ite-its+1 ) :: sols ! Downward solar rad onto surface (sw direct)
  359. real(r8), dimension( 1:ite-its+1 ) :: soll ! Downward solar rad onto surface (lw direct)
  360. real(r8), dimension( 1:ite-its+1 ) :: solsd ! Downward solar rad onto surface (sw diffuse)
  361. real(r8), dimension( 1:ite-its+1 ) :: solld ! Downward solar rad onto surface (lw diffuse)
  362. real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: qrs ! Solar heating rate
  363. real(r8), dimension( 1:ite-its+1 ) :: fsds ! Flux Shortwave Downwelling Surface
  364. real(r8), dimension( 1:ite-its+1, 1:kte-kts+1 ) :: qrl ! Longwave cooling rate
  365. real(r8), dimension( 1:ite-its+1 ) :: flwds ! Surface down longwave flux
  366. real(r8), dimension( 1:ite-its+1, levsiz, num_months ) :: ozmixmj ! monthly ozone mixing ratio
  367. real(r8), dimension( 1:ite-its+1, levsiz ) :: ozmix ! ozone mixing ratio (time interpolated)
  368. real(r8), dimension(levsiz) :: pin ! ozone pressure level
  369. real(r8), dimension(1:ite-its+1) :: m_psjp,m_psjn ! MATCH surface pressure
  370. real(r8), dimension( 1:ite-its+1, paerlev, naer_c ) :: aerosoljp ! monthly aerosol concentrations
  371. real(r8), dimension( 1:ite-its+1, paerlev, naer_c ) :: aerosoljn ! monthly aerosol concentrations
  372. real(r8), dimension(paerlev) :: m_hybi
  373. real(r8), dimension(1:ite-its+1 ) :: clat ! latitude in radians for columns
  374. real(r8), dimension(its:ite,kts:kte+1,kts:kte+1) :: abstot ! Total absorptivity
  375. real(r8), dimension(its:ite,kts:kte,4) :: absnxt ! Total nearest layer absorptivity
  376. real(r8), dimension(its:ite,kts:kte+1) :: emstot ! Total emissivity
  377. CHARACTER(LEN=256) :: msgstr
  378. #if !defined(MAC_KLUDGE)
  379. lchnk = 1
  380. begchunk = ims
  381. endchunk = ime
  382. ncol = ite - its + 1
  383. pcols= ite - its + 1
  384. pver = kte - kts + 1
  385. pverp= pver + 1
  386. pverr = kte - kts + 1
  387. pverrp= pverr + 1
  388. ! number of advected constituents and non-advected constituents (including water vapor)
  389. ppcnst = n_cldadv
  390. ! number of non-advected constituents
  391. pnats = 0
  392. pcnst = ppcnst-pnats
  393. ! check the # species defined for the input climatology and naer
  394. ! if(naer_c.ne.naer) then
  395. ! WRITE( wrf_err_message , * ) 'naer_c ne naer ', naer_c, naer
  396. if(naer_c.ne.naer_all) then
  397. WRITE( wrf_err_message , * ) 'naer_c-1 ne naer_all ', naer_c, naer_all
  398. CALL wrf_error_fatal ( wrf_err_message )
  399. endif
  400. ! update CO2 volume mixing ratio (co2vmr)
  401. ! determine time interpolation factors, check sanity
  402. ! of interpolation factors to within 32-bit roundoff
  403. ! assume that day of year is 1 for all input data
  404. !
  405. nyrm = yr - yrdata(1) + 1
  406. nyrp = nyrm + 1
  407. doymodel = yr*365. + julian
  408. doydatam = yrdata(nyrm)*365. + 1.
  409. doydatap = yrdata(nyrp)*365. + 1.
  410. deltat = doydatap - doydatam
  411. fact1 = (doydatap - doymodel)/deltat
  412. fact2 = (doymodel - doydatam)/deltat
  413. co2vmr = (co2(nyrm)*fact1 + co2(nyrp)*fact2)*1.e-06
  414. co2mmr=co2vmr*mwco2/mwdry
  415. !
  416. !===================================================
  417. ! Radiation computations
  418. !===================================================
  419. do k=1,levsiz
  420. pin(k)=pin0(k)
  421. enddo
  422. do k=1,paerlev
  423. m_hybi(k)=m_hybi0(k)
  424. enddo
  425. ! check for uninitialized arrays
  426. if(abstot_3d(its,kts,kts,jts) .eq. 0.0 .and. .not.doabsems .and. dolw)then
  427. CALL wrf_debug(0, 'camrad lw: CAUTION: re-calculating abstot, absnxt, emstot on restart')
  428. doabsems = .true.
  429. endif
  430. do j =jts,jte
  431. !
  432. ! Cosine solar zenith angle for current time step
  433. !
  434. ! call zenith (calday, clat, clon, coszrs, ncol)
  435. do i = its,ite
  436. ii = i - its + 1
  437. ! XT24 is the fractional part of simulation days plus half of RADT expressed in
  438. ! units of minutes
  439. ! JULIAN is in days
  440. ! RADT is in minutes
  441. XT24=MOD(XTIME+RADT*0.5,1440.)
  442. TLOCTM=GMT+XT24/60.+XLONG(I,J)/15.
  443. HRANG=15.*(TLOCTM-12.)*DEGRAD
  444. XXLAT=XLAT(I,J)*DEGRAD
  445. clat(ii)=xxlat
  446. coszrs(II)=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG)
  447. enddo
  448. ! moist variables
  449. do k = kts,kte
  450. kk = kte - k + kts
  451. do i = its,ite
  452. ii = i - its + 1
  453. ! convert to specific humidity
  454. q(ii,kk,1) = max(1.e-10,qv3d(i,k,j)/(1.+qv3d(i,k,j)))
  455. IF ( F_QI .and. F_QC .and. F_QS ) THEN
  456. q(ii,kk,ixcldliq) = max(0.,qc3d(i,k,j)/(1.+qv3d(i,k,j)))
  457. q(ii,kk,ixcldice) = max(0.,(qi3d(i,k,j)+qs3d(i,k,j))/(1.+qv3d(i,k,j)))
  458. ELSE IF ( F_QC .and. F_QR ) THEN
  459. ! Warm rain or simple ice
  460. q(ii,kk,ixcldliq) = 0.
  461. q(ii,kk,ixcldice) = 0.
  462. if(t_phy(i,k,j).gt.273.15)q(ii,kk,ixcldliq) = max(0.,qc3d(i,k,j)/(1.+qv3d(i,k,j)))
  463. if(t_phy(i,k,j).le.273.15)q(ii,kk,ixcldice) = max(0.,qc3d(i,k,j)/(1.+qv3d(i,k,j)))
  464. ELSE IF ( F_QC .and. F_QS ) THEN
  465. ! For Ferrier (note that currently Ferrier has QI, so this section will not be used)
  466. q(ii,kk,ixcldice) = max(0.,qc3d(i,k,j)/(1.+qv3d(i,k,j))*f_ice_phy(i,k,j))
  467. 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)))
  468. ELSE
  469. q(ii,kk,ixcldliq) = 0.
  470. q(ii,kk,ixcldice) = 0.
  471. ENDIF
  472. cld(ii,kk) = CLDFRA(I,K,J)
  473. enddo
  474. enddo
  475. do i = its,ite
  476. ii = i - its + 1
  477. landfrac(ii) = 2.-XLAND(I,J)
  478. landm(ii) = landfrac(ii)
  479. snowh(ii) = 0.001*SNOW(I,J)
  480. icefrac(ii) = XICE(I,J)
  481. enddo
  482. do m=1,num_months-1
  483. do k=1,levsiz
  484. do i = its,ite
  485. ii = i - its + 1
  486. ozmixmj(ii,k,m) = ozmixm(i,k,j,m+1)
  487. enddo
  488. enddo
  489. enddo
  490. do i = its,ite
  491. ii = i - its + 1
  492. m_psjp(ii) = m_psp(i,j)
  493. m_psjn(ii) = m_psn(i,j)
  494. enddo
  495. do n=1,naer_c
  496. do k=1,paerlev
  497. do i = its,ite
  498. ii = i - its + 1
  499. aerosoljp(ii,k,n) = aerosolcp(i,k,j,n)
  500. aerosoljn(ii,k,n) = aerosolcn(i,k,j,n)
  501. enddo
  502. enddo
  503. enddo
  504. !
  505. ! Complete radiation calculations
  506. !
  507. do i = its,ite
  508. ii = i - its + 1
  509. lwups(ii) = stebol*EMISS(I,J)*TSK(I,J)**4
  510. enddo
  511. do k = kts,kte+1
  512. kk = kte - k + kts + 1
  513. do i = its,ite
  514. ii = i - its + 1
  515. pint(ii,kk) = p8w(i,k,j)
  516. if(k.eq.kts)ps(ii)=pint(ii,kk)
  517. lnpint(ii,kk) = log(pint(ii,kk))
  518. enddo
  519. enddo
  520. if(.not.doabsems .and. dolw)then
  521. ! do kk = kts,kte+1
  522. do kk = 1,cam_abs_dim2
  523. do kk1 = kts,kte+1
  524. do i = its,ite
  525. abstot(i,kk1,kk) = abstot_3d(i,kk1,kk,j)
  526. enddo
  527. enddo
  528. enddo
  529. ! do kk = 1,4
  530. do kk = 1,cam_abs_dim1
  531. do kk1 = kts,kte
  532. do i = its,ite
  533. absnxt(i,kk1,kk) = absnxt_3d(i,kk1,kk,j)
  534. enddo
  535. enddo
  536. enddo
  537. do kk = kts,kte+1
  538. do i = its,ite
  539. emstot(i,kk) = emstot_3d(i,kk,j)
  540. enddo
  541. enddo
  542. endif
  543. do k = kts,kte
  544. kk = kte - k + kts
  545. do i = its,ite
  546. ii = i - its + 1
  547. pmid(ii,kk) = p_phy(i,k,j)
  548. lnpmid(ii,kk) = log(pmid(ii,kk))
  549. lnpint(ii,kk) = log(pint(ii,kk))
  550. pdel(ii,kk) = pint(ii,kk+1) - pint(ii,kk)
  551. t(ii,kk) = t_phy(i,k,j)
  552. zm(ii,kk) = z(i,k,j)
  553. enddo
  554. enddo
  555. ! Compute cloud water/ice paths and optical properties for input to radiation
  556. call param_cldoptics_calc(ncol, pcols, pver, pverp, pverr, pverrp, ppcnst, q, cld, landfrac, landm,icefrac, &
  557. pdel, t, ps, pmid, pint, cicewp, cliqwp, emis, rel, rei, pmxrgn, nmxrgn, snowh)
  558. !-----fds (06/2010)----------------------------
  559. SELECT CASE(sf_surface_physics)
  560. CASE (SSIBSCHEME)
  561. if (xtime .gt. 1.0) then
  562. ! call wrf_message("using SSiB albedoes for land points")
  563. do i = its,ite
  564. ii = i - its + 1
  565. if (xland(i,j).lt.1.5) then !land points only
  566. asdir(ii) = ALSWVISDIR(i,j) ! SSiB visdir albedo
  567. asdif(ii) = ALSWVISDIF(i,j) ! SSiB visdif albedo
  568. aldir(ii) = ALSWNIRDIR(i,j) ! SSiB nirdir albedo
  569. aldif(ii) = ALSWNIRDIF(i,j) ! SSiB nirdif albedo
  570. else
  571. asdir(ii) = albedo(i,j)
  572. asdif(ii) = albedo(i,j)
  573. aldir(ii) = albedo(i,j)
  574. aldif(ii) = albedo(i,j)
  575. endif
  576. enddo
  577. else
  578. do i = its,ite
  579. ii = i - its + 1
  580. asdir(ii) = albedo(i,j)
  581. asdif(ii) = albedo(i,j)
  582. aldir(ii) = albedo(i,j)
  583. aldif(ii) = albedo(i,j)
  584. enddo
  585. endif
  586. CASE DEFAULT
  587. do i = its,ite
  588. ii = i - its + 1
  589. ! use same albedo for direct and diffuse
  590. ! change this when separate values are provided
  591. asdir(ii) = albedo(i,j)
  592. asdif(ii) = albedo(i,j)
  593. aldir(ii) = albedo(i,j)
  594. aldif(ii) = albedo(i,j)
  595. enddo
  596. END SELECT
  597. !-----------------------------------------------
  598. ! WRF allocate space here (not needed if oznini is called)
  599. ! allocate (ozmix(pcols,levsiz,begchunk:endchunk)) ! This line from oznini.F90
  600. call radctl (j,lchnk, ncol, pcols, pver, pverp, pverr, pverrp, ppcnst, pcnst, lwups, emis, pmid, &
  601. pint, lnpmid, lnpint, pdel, t, q, &
  602. cld, cicewp, cliqwp, tauxcl, tauxci, coszrs, clat, asdir, asdif, &
  603. aldir, aldif, solcon, GMT,JULDAY,JULIAN,DT,XTIME, &
  604. pin, ozmixmj, ozmix, levsiz, num_months, &
  605. m_psjp,m_psjn, aerosoljp, aerosoljn, m_hybi, paerlev, naer_c, pmxrgn, nmxrgn, &
  606. dolw, dosw, doabsems, abstot, absnxt, emstot, &
  607. fsup, fsupc, fsdn, fsdnc, flup, flupc, fldn, fldnc, swcftoa, lwcftoa, olrtoa, &
  608. fsns, fsnt ,flns ,flnt , &
  609. qrs, qrl, flwds, rel, rei, &
  610. sols, soll, solsd, solld, &
  611. landfrac, zm, fsds)
  612. do k = kts,kte
  613. kk = kte - k + kts
  614. do i = its,ite
  615. ii = i - its + 1
  616. if(dolw)RTHRATENLW(I,K,J) = 1.e4*qrl(ii,kk)/(cpair*pi_phy(i,k,j))
  617. if(dosw)RTHRATENSW(I,K,J) = 1.e4*qrs(ii,kk)/(cpair*pi_phy(i,k,j))
  618. cemiss(i,k,j) = emis(ii,kk)
  619. taucldc(i,k,j) = tauxcl(ii,kk)
  620. taucldi(i,k,j) = tauxci(ii,kk)
  621. enddo
  622. enddo
  623. if(doabsems .and. dolw)then
  624. ! do kk = kts,kte+1
  625. do kk = 1,cam_abs_dim2
  626. do kk1 = kts,kte+1
  627. do i = its,ite
  628. abstot_3d(i,kk1,kk,j) = abstot(i,kk1,kk)
  629. enddo
  630. enddo
  631. enddo
  632. ! do kk = 1,4
  633. do kk = 1,cam_abs_dim1
  634. do kk1 = kts,kte
  635. do i = its,ite
  636. absnxt_3d(i,kk1,kk,j) = absnxt(i,kk1,kk)
  637. enddo
  638. enddo
  639. enddo
  640. do kk = kts,kte+1
  641. do i = its,ite
  642. emstot_3d(i,kk,j) = emstot(i,kk)
  643. enddo
  644. enddo
  645. endif
  646. IF(PRESENT(SWUPT))THEN
  647. if(dosw)then
  648. ! Added shortwave and longwave upward/downward total and clear sky fluxes
  649. do k = kts,kte+1
  650. kk = kte +1 - k + kts
  651. do i = its,ite
  652. ii = i - its + 1
  653. ! swup(i,k,j) = fsup(ii,kk)
  654. ! swupclear(i,k,j) = fsupc(ii,kk)
  655. ! swdn(i,k,j) = fsdn(ii,kk)
  656. ! swdnclear(i,k,j) = fsdnc(ii,kk)
  657. if(k.eq.kte+1)then
  658. swupt(i,j) = fsup(ii,kk)
  659. swuptc(i,j) = fsupc(ii,kk)
  660. swdnt(i,j) = fsdn(ii,kk)
  661. swdntc(i,j) = fsdnc(ii,kk)
  662. endif
  663. if(k.eq.kts)then
  664. swupb(i,j) = fsup(ii,kk)
  665. swupbc(i,j) = fsupc(ii,kk)
  666. swdnb(i,j) = fsdn(ii,kk)
  667. swdnbc(i,j) = fsdnc(ii,kk)
  668. endif
  669. ! if(i.eq.30.and.j.eq.30) then
  670. ! print 1234, 'short ', i,ii,k,kk,fsup(ii,kk),fsupc(ii,kk),fsdn(ii,kk),fsdnc(ii,kk)
  671. ! 1234 format (a6,4i4,4f10.3)
  672. ! endif
  673. enddo
  674. enddo
  675. endif
  676. if(dolw)then
  677. ! Added shortwave and longwave upward/downward total and clear sky fluxes
  678. do k = kts,kte+1
  679. kk = kte +1 - k + kts
  680. do i = its,ite
  681. ii = i - its + 1
  682. ! lwup(i,k,j) = flup(ii,kk)
  683. ! lwupclear(i,k,j) = flupc(ii,kk)
  684. ! lwdn(i,k,j) = fldn(ii,kk)
  685. ! lwdnclear(i,k,j) = fldnc(ii,kk)
  686. if(k.eq.kte+1)then
  687. lwupt(i,j) = flup(ii,kk)
  688. lwuptc(i,j) = flupc(ii,kk)
  689. lwdnt(i,j) = fldn(ii,kk)
  690. lwdntc(i,j) = fldnc(ii,kk)
  691. endif
  692. if(k.eq.kts)then
  693. lwupb(i,j) = flup(ii,kk)
  694. lwupbc(i,j) = flupc(ii,kk)
  695. lwdnb(i,j) = fldn(ii,kk)
  696. lwdnbc(i,j) = fldnc(ii,kk)
  697. endif
  698. ! if(i.eq.30.and.j.eq.30) then
  699. ! print 1234, 'long ', i,ii,k,kk,flup(ii,kk),flupc(ii,kk),fldn(ii,kk),fldnc(ii,kk)
  700. ! 1234 format (a6,4i4,4f10.3)
  701. ! endif
  702. enddo
  703. enddo
  704. endif
  705. ENDIF
  706. do i = its,ite
  707. ii = i - its + 1
  708. ! Added shortwave and longwave cloud forcing at TOA and surface
  709. if(dolw)then
  710. GLW(I,J) = flwds(ii)
  711. lwcf(i,j) = lwcftoa(ii)
  712. olr(i,j) = olrtoa(ii)
  713. endif
  714. if(dosw)then
  715. GSW(I,J) = fsns(ii)
  716. swcf(i,j) = swcftoa(ii)
  717. coszr(i,j) = coszrs(ii)
  718. endif
  719. enddo
  720. !-------fds (06/2010)---------
  721. SELECT CASE(sf_surface_physics)
  722. CASE (SSIBSCHEME)
  723. ! call wrf_message("CAM using ssib albedo2")
  724. if(dosw)then
  725. do i = its,ite
  726. ii = i - its + 1
  727. SWVISDIR(I,J) = sols(ii) !SSiB
  728. SWVISDIF(I,J) = solsd(ii) !SSiB
  729. SWNIRDIR(I,J) = soll(ii) !SSiB
  730. SWNIRDIF(I,J) = solld(ii) !SSiB
  731. enddo
  732. endif
  733. END SELECT
  734. !-----------------------------
  735. enddo ! j-loop
  736. #endif
  737. end subroutine camrad
  738. !====================================================================
  739. SUBROUTINE camradinit( &
  740. R_D,R_V,CP,G,STBOLT,EP_2,shalf,pptop, &
  741. ozmixm,pin,levsiz,XLAT,num_months, &
  742. m_psp,m_psn,m_hybi,aerosolcp,aerosolcn, &
  743. paerlev,naer_c, &
  744. ids, ide, jds, jde, kds, kde, &
  745. ims, ime, jms, jme, kms, kme, &
  746. its, ite, jts, jte, kts, kte )
  747. USE module_wrf_error
  748. USE module_state_description
  749. !USE module_configure
  750. !--------------------------------------------------------------------
  751. IMPLICIT NONE
  752. !--------------------------------------------------------------------
  753. INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
  754. ims, ime, jms, jme, kms, kme, &
  755. its, ite, jts, jte, kts, kte
  756. REAL, intent(in) :: pptop
  757. REAL, INTENT(IN) :: R_D,R_V,CP,G,STBOLT,EP_2
  758. REAL, DIMENSION( kms:kme ) :: shalf
  759. INTEGER, INTENT(IN ) :: levsiz, num_months
  760. INTEGER, INTENT(IN ) :: paerlev, naer_c
  761. REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: XLAT
  762. REAL, DIMENSION( ims:ime, levsiz, jms:jme, num_months ), &
  763. INTENT(INOUT ) :: OZMIXM
  764. REAL, DIMENSION(levsiz), INTENT(INOUT ) :: PIN
  765. REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT ) :: m_psp,m_psn
  766. REAL, DIMENSION(paerlev), INTENT(INOUT ) :: m_hybi
  767. REAL, DIMENSION( ims:ime, paerlev, jms:jme, naer_c ), &
  768. INTENT(INOUT) :: aerosolcp,aerosolcn
  769. REAL(r8) :: pstd
  770. REAL(r8) :: rh2o, cpair
  771. ! These were made allocatable 20090612 to save static memory allocation. JM
  772. IF ( .NOT. ALLOCATED( ksul ) ) ALLOCATE( ksul( nrh, nspint ) )
  773. IF ( .NOT. ALLOCATED( wsul ) ) ALLOCATE( wsul( nrh, nspint ) )
  774. IF ( .NOT. ALLOCATED( gsul ) ) ALLOCATE( gsul( nrh, nspint ) )
  775. IF ( .NOT. ALLOCATED( ksslt ) ) ALLOCATE( ksslt( nrh, nspint ) )
  776. IF ( .NOT. ALLOCATED( wsslt ) ) ALLOCATE( wsslt( nrh, nspint ) )
  777. IF ( .NOT. ALLOCATED( gsslt ) ) ALLOCATE( gsslt( nrh, nspint ) )
  778. IF ( .NOT. ALLOCATED( kcphil ) ) ALLOCATE( kcphil( nrh, nspint ) )
  779. IF ( .NOT. ALLOCATED( wcphil ) ) ALLOCATE( wcphil( nrh, nspint ) )
  780. IF ( .NOT. ALLOCATED( gcphil ) ) ALLOCATE( gcphil( nrh, nspint ) )
  781. IF ( .NOT. ALLOCATED(ah2onw ) ) ALLOCATE( ah2onw(n_p, n_tp, n_u, n_te, n_rh) )
  782. IF ( .NOT. ALLOCATED(eh2onw ) ) ALLOCATE( eh2onw(n_p, n_tp, n_u, n_te, n_rh) )
  783. IF ( .NOT. ALLOCATED(ah2ow ) ) ALLOCATE( ah2ow(n_p, n_tp, n_u, n_te, n_rh) )
  784. IF ( .NOT. ALLOCATED(cn_ah2ow) ) ALLOCATE( cn_ah2ow(n_p, n_tp, n_u, n_te, n_rh) )
  785. IF ( .NOT. ALLOCATED(cn_eh2ow) ) ALLOCATE( cn_eh2ow(n_p, n_tp, n_u, n_te, n_rh) )
  786. IF ( .NOT. ALLOCATED(ln_ah2ow) ) ALLOCATE( ln_ah2ow(n_p, n_tp, n_u, n_te, n_rh) )
  787. IF ( .NOT. ALLOCATED(ln_eh2ow) ) ALLOCATE( ln_eh2ow(n_p, n_tp, n_u, n_te, n_rh) )
  788. #if !defined(MAC_KLUDGE)
  789. ozncyc = .true.
  790. indirect = .true.
  791. ixcldliq = 2
  792. ixcldice = 3
  793. #if (NMM_CORE != 1)
  794. ! aerosol array is not in the NMM Registry
  795. ! since CAM radiation not available to NMM (yet)
  796. ! so this is blocked out to enable CAM compilation with NMM
  797. idxSUL = P_SUL
  798. idxSSLT = P_SSLT
  799. idxDUSTfirst = P_DUST1
  800. idxOCPHO = P_OCPHO
  801. idxCARBONfirst = P_OCPHO
  802. idxBCPHO = P_BCPHO
  803. idxOCPHI = P_OCPHI
  804. idxBCPHI = P_BCPHI
  805. idxBG = P_BG
  806. idxVOLC = P_VOLC
  807. #endif
  808. pstd = 101325.0
  809. ! from physconst module
  810. mwdry = 28.966 ! molecular weight dry air ~ kg/kmole (shr_const_mwdair)
  811. mwco2 = 44. ! molecular weight co2
  812. mwh2o = 18.016 ! molecular weight water vapor (shr_const_mwwv)
  813. mwch4 = 16. ! molecular weight ch4
  814. mwn2o = 44. ! molecular weight n2o
  815. mwf11 = 136. ! molecular weight cfc11
  816. mwf12 = 120. ! molecular weight cfc12
  817. cappa = R_D/CP
  818. rair = R_D
  819. tmelt = 273.16 ! freezing T of fresh water ~ K
  820. r_universal = 6.02214e26 * STBOLT ! Universal gas constant ~ J/K/kmole
  821. latvap = 2.501e6 ! latent heat of evaporation ~ J/kg
  822. latice = 3.336e5 ! latent heat of fusion ~ J/kg
  823. zvir = R_V/R_D - 1.
  824. rh2o = R_V
  825. cpair = CP
  826. !
  827. epsqs = EP_2
  828. CALL radini(G, CP, EP_2, STBOLT, pstd*10.0 )
  829. CALL esinti(epsqs ,latvap ,latice ,rh2o ,cpair ,tmelt )
  830. CALL oznini(ozmixm,pin,levsiz,num_months,XLAT, &
  831. ids, ide, jds, jde, kds, kde, &
  832. ims, ime, jms, jme, kms, kme, &
  833. its, ite, jts, jte, kts, kte)
  834. CALL aerosol_init(m_psp,m_psn,m_hybi,aerosolcp,aerosolcn,paerlev,naer_c,shalf,pptop, &
  835. ids, ide, jds, jde, kds, kde, &
  836. ims, ime, jms, jme, kms, kme, &
  837. its, ite, jts, jte, kts, kte)
  838. #endif
  839. END SUBROUTINE camradinit
  840. #if !defined(MAC_KLUDGE)
  841. subroutine oznint(julday,julian,dt,gmt,xtime,ozmixmj,ozmix,levsiz,num_months,pcols)
  842. IMPLICIT NONE
  843. INTEGER, INTENT(IN ) :: levsiz, num_months,pcols
  844. REAL(r8), DIMENSION( pcols, levsiz, num_months ), &
  845. INTENT(IN ) :: ozmixmj
  846. REAL, INTENT(IN ) :: XTIME,GMT
  847. INTEGER, INTENT(IN ) :: JULDAY
  848. REAL, INTENT(IN ) :: JULIAN
  849. REAL, INTENT(IN ) :: DT
  850. REAL(r8), DIMENSION( pcols, levsiz ), &
  851. INTENT(OUT ) :: ozmix
  852. !Local
  853. REAL(r8) :: intJULIAN
  854. integer :: np1,np,nm,m,k,i
  855. integer :: IJUL
  856. integer, dimension(12) :: date_oz
  857. data date_oz/16, 45, 75, 105, 136, 166, 197, 228, 258, 289, 319, 350/
  858. real(r8) :: cdayozp, cdayozm
  859. real(r8) :: fact1, fact2
  860. logical :: finddate
  861. CHARACTER(LEN=256) :: msgstr
  862. ! JULIAN starts from 0.0 at 0Z on 1 Jan.
  863. intJULIAN = JULIAN + 1.0_r8 ! offset by one day
  864. ! jan 1st 00z is julian=1.0 here
  865. IJUL=INT(intJULIAN)
  866. ! Note that following will drift.
  867. ! Need to use actual month/day info to compute julian.
  868. intJULIAN=intJULIAN-FLOAT(IJUL)
  869. IJUL=MOD(IJUL,365)
  870. IF(IJUL.EQ.0)IJUL=365
  871. intJULIAN=intJULIAN+IJUL
  872. np1=1
  873. finddate=.false.
  874. ! do m=1,num_months
  875. do m=1,12
  876. if(date_oz(m).gt.intjulian.and..not.finddate) then
  877. np1=m
  878. finddate=.true.
  879. endif
  880. enddo
  881. cdayozp=date_oz(np1)
  882. if(np1.gt.1) then
  883. cdayozm=date_oz(np1-1)
  884. np=np1
  885. nm=np-1
  886. else
  887. cdayozm=date_oz(12)
  888. np=np1
  889. nm=12
  890. endif
  891. call getfactors(ozncyc,np1, cdayozm, cdayozp,intjulian, &
  892. fact1, fact2)
  893. !
  894. ! Time interpolation.
  895. !
  896. do k=1,levsiz
  897. do i=1,pcols
  898. ozmix(i,k) = ozmixmj(i,k,nm)*fact1 + ozmixmj(i,k,np)*fact2
  899. end do
  900. end do
  901. END subroutine oznint
  902. subroutine get_aerosol(c, julday, julian, dt, gmt, xtime, m_psp, m_psn, aerosoljp, &
  903. aerosoljn, m_hybi, paerlev, naer_c, pint, pcols, pver, pverp, pverr, pverrp, AEROSOLt, scale)
  904. !------------------------------------------------------------------
  905. !
  906. ! Input:
  907. ! time at which aerosol mmrs are needed (get_curr_calday())
  908. ! chunk index
  909. ! CAM's vertical grid (pint)
  910. !
  911. ! Output:
  912. ! values for Aerosol Mass Mixing Ratios at specified time
  913. ! on vertical grid specified by CAM (AEROSOLt)
  914. !
  915. ! Method:
  916. ! first determine which indexs of aerosols are the bounding data sets
  917. ! interpolate both onto vertical grid aerm(),aerp().
  918. ! from those two, interpolate in time.
  919. !
  920. !------------------------------------------------------------------
  921. ! use volcanicmass, only: get_volcanic_mass
  922. ! use timeinterp, only: getfactors
  923. !
  924. ! aerosol fields interpolated to current time step
  925. ! on pressure levels of this time step.
  926. ! these should be made read-only for other modules
  927. ! Is allocation done correctly here?
  928. !
  929. integer, intent(in) :: c ! Chunk Id.
  930. integer, intent(in) :: paerlev, naer_c, pcols, pver, pverp, pverr, pverrp
  931. real(r8), intent(in) :: pint(pcols,pverp) ! midpoint pres.
  932. real(r8), intent(in) :: scale(naer_all) ! scale each aerosol by this amount
  933. REAL, INTENT(IN ) :: XTIME,GMT
  934. INTEGER, INTENT(IN ) :: JULDAY
  935. REAL, INTENT(IN ) :: JULIAN
  936. REAL, INTENT(IN ) :: DT
  937. real(r8), intent(in ) :: m_psp(pcols),m_psn(pcols) ! Match surface pressure
  938. real(r8), intent(in ) :: aerosoljp(pcols,paerlev,naer_c)
  939. real(r8), intent(in ) :: aerosoljn(pcols,paerlev,naer_c)
  940. real(r8), intent(in ) :: m_hybi(paerlev)
  941. real(r8), intent(out) :: AEROSOLt(pcols, pver, naer_all) ! aerosols
  942. !
  943. ! Local workspace
  944. !
  945. real(r8) caldayloc ! calendar day of current timestep
  946. real(r8) fact1, fact2 ! time interpolation factors
  947. integer :: nm = 1 ! index to prv month in array. init to 1 and toggle between 1 and 2
  948. integer :: np = 2 ! index to nxt month in array. init to 2 and toggle between 1 and 2
  949. integer :: mo_nxt = bigint ! index to nxt month in file
  950. integer :: mo_prv ! index to previous month
  951. real(r8) :: cdaym = inf ! calendar day of prv month
  952. real(r8) :: cdayp = inf ! calendar day of next month
  953. real(r8) :: Mid(12) ! Days into year for mid month date
  954. 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 /
  955. integer i, k, j ! spatial indices
  956. integer m ! constituent index
  957. integer lats(pcols),lons(pcols) ! latitude and longitudes of column
  958. integer ncol ! number of columns
  959. INTEGER IJUL
  960. REAL(r8) intJULIAN
  961. real(r8) speciesmin(naer) ! minimal value for each species
  962. !
  963. ! values before current time step "the minus month"
  964. ! aerosolm(pcols,pver) is value of preceeding month's aerosol mmr
  965. ! aerosolp(pcols,pver) is value of next month's aerosol mmr
  966. ! (think minus and plus or values to left and right of point to be interpolated)
  967. !
  968. real(r8) AEROSOLm(pcols,pver,naer) ! aerosol mmr from MATCH in column at previous (minus) month
  969. !
  970. ! values beyond (or at) current time step "the plus month"
  971. !
  972. real(r8) AEROSOLp(pcols,pver,naer) ! aerosol mmr from MATCH in column at next (plus) month
  973. CHARACTER(LEN=256) :: msgstr
  974. ! JULIAN starts from 0.0 at 0Z on 1 Jan.
  975. intJULIAN = JULIAN + 1.0_r8 ! offset by one day
  976. ! jan 1st 00z is julian=1.0 here
  977. IJUL=INT(intJULIAN)
  978. ! Note that following will drift.
  979. ! Need to use actual month/day info to compute julian.
  980. intJULIAN=intJULIAN-FLOAT(IJUL)
  981. IJUL=MOD(IJUL,365)
  982. IF(IJUL.EQ.0)IJUL=365
  983. caldayloc=intJULIAN+IJUL
  984. if (caldayloc < Mid(1)) then
  985. mo_prv = 12
  986. mo_nxt = 1
  987. else if (caldayloc >= Mid(12)) then
  988. mo_prv = 12
  989. mo_nxt = 1
  990. else
  991. do i = 2 , 12
  992. if (caldayloc < Mid(i)) then
  993. mo_prv = i-1
  994. mo_nxt = i
  995. exit
  996. end if
  997. end do
  998. end if
  999. !
  1000. ! Set initial calendar day values
  1001. !
  1002. cdaym = Mid(mo_prv)
  1003. cdayp = Mid(mo_nxt)
  1004. !
  1005. ! Determine time interpolation factors. 1st arg says we are cycling 1 year of data
  1006. !
  1007. call getfactors (.true., mo_nxt, cdaym, cdayp, caldayloc, &
  1008. fact1, fact2)
  1009. !
  1010. ! interpolate (prv and nxt month) bounding datasets onto cam vertical grid.
  1011. ! compute mass mixing ratios on CAMS's pressure coordinate
  1012. ! for both the "minus" and "plus" months
  1013. !
  1014. ! ncol = get_ncols_p(c)
  1015. ncol = pcols
  1016. ! call vert_interpolate (M_ps_cam_col(1,c,nm), pint, nm, AEROSOLm, ncol, c)
  1017. ! call vert_interpolate (M_ps_cam_col(1,c,np), pint, np, AEROSOLp, ncol, c)
  1018. call vert_interpolate (m_psp, aerosoljp, m_hybi, paerlev, naer_c, pint, nm, AEROSOLm, pcols, pver, pverp, ncol, c)
  1019. call vert_interpolate (m_psn, aerosoljn, m_hybi, paerlev, naer_c, pint, np, AEROSOLp, pcols, pver, pverp, ncol, c)
  1020. !
  1021. ! Time interpolate.
  1022. !
  1023. do m=1,naer
  1024. do k=1,pver
  1025. do i=1,ncol
  1026. AEROSOLt(i,k,m) = AEROSOLm(i,k,m)*fact1 + AEROSOLp(i,k,m)*fact2
  1027. end do
  1028. end do
  1029. end do
  1030. ! do i=1,ncol
  1031. ! Match_ps_chunk(i,c) = m_ps(i,nm)*fact1 + m_ps(i,np)*fact2
  1032. ! end do
  1033. !
  1034. ! get background aerosol (tuning) field
  1035. !
  1036. call background (c, ncol, pint, pcols, pverr, pverrp, AEROSOLt(:, :, idxBG))
  1037. !
  1038. ! find volcanic aerosol masses
  1039. !
  1040. ! if (strat_volcanic) then
  1041. ! call get_volcanic_mass(c, AEROSOLt(:,:,idxVOLC))
  1042. ! else
  1043. AEROSOLt(:,:,idxVOLC) = 0._r8
  1044. ! endif
  1045. !
  1046. ! exit if mmr is negative (we have previously set
  1047. ! cumulative mass to be a decreasing function.)
  1048. !
  1049. speciesmin(:) = 0. ! speciesmin(m) = 0 is minimum mmr for each species
  1050. do m=1,naer
  1051. do k=1,pver
  1052. do i=1,ncol
  1053. if (AEROSOLt(i, k, m) < speciesmin(m)) then
  1054. write(6,*) 'AEROSOL_INTERPOLATE: negative mass mixing ratio, exiting'
  1055. write(6,*) 'm, column, pver',m, i, k ,AEROSOLt(i, k, m)
  1056. call endrun ()
  1057. end if
  1058. end do
  1059. end do
  1060. end do
  1061. !
  1062. ! scale any AEROSOLS as required
  1063. !
  1064. call scale_aerosols (AEROSOLt, pcols, pver, ncol, c, scale)
  1065. return
  1066. end subroutine get_aerosol
  1067. subroutine aerosol_indirect(ncol,lchnk,pcols,pver,ppcnst,landfrac,pmid,t,qm1,cld,zm,rel)
  1068. !--------------------------------------------------------------
  1069. ! Compute effect of sulfate on effective liquid water radius
  1070. ! Method of Martin et. al.
  1071. !--------------------------------------------------------------
  1072. ! use constituents, only: ppcnst, cnst_get_ind
  1073. ! use history, only: outfld
  1074. !#include <comctl.h>
  1075. integer, intent(in) :: ncol ! number of atmospheric columns
  1076. integer, intent(in) :: lchnk ! chunk identifier
  1077. integer, intent(in) :: pcols,pver,ppcnst
  1078. real(r8), intent(in) :: landfrac(pcols) ! land fraction
  1079. real(r8), intent(in) :: pmid(pcols,pver) ! Model level pressures
  1080. real(r8), intent(in) :: t(pcols,pver) ! Model level temperatures
  1081. real(r8), intent(in) :: qm1(pcols,pver,ppcnst) ! Specific humidity and tracers
  1082. real(r8), intent(in) :: cld(pcols,pver) ! Fractional cloud cover
  1083. real(r8), intent(in) :: zm(pcols,pver) ! Height of midpoints (above surface)
  1084. real(r8), intent(in) :: rel(pcols,pver) ! liquid effective drop size (microns)
  1085. !
  1086. ! local variables
  1087. !
  1088. real(r8) locrhoair(pcols,pver) ! dry air density [kg/m^3 ]
  1089. real(r8) lwcwat(pcols,pver) ! in-cloud liquid water path [kg/m^3 ]
  1090. real(r8) sulfmix(pcols,pver) ! sulfate mass mixing ratio [kg/kg ]
  1091. real(r8) so4mass(pcols,pver) ! sulfate mass concentration [g/cm^3 ]
  1092. real(r8) Aso4(pcols,pver) ! sulfate # concentration [#/cm^3 ]
  1093. real(r8) Ntot(pcols,pver) ! ccn # concentration [#/cm^3 ]
  1094. real(r8) relmod(pcols,pver) ! effective radius [microns]
  1095. real(r8) wrel(pcols,pver) ! weighted effective radius [microns]
  1096. real(r8) wlwc(pcols,pver) ! weighted liq. water content [kg/m^3 ]
  1097. real(r8) cldfrq(pcols,pver) ! frequency of occurance of...
  1098. ! ! clouds (cld => 0.01) [fraction]
  1099. real(r8) locPi ! my piece of the pi
  1100. real(r8) Rdryair ! gas constant of dry air [J/deg/kg]
  1101. real(r8) rhowat ! density of water [kg/m^3 ]
  1102. re

Large files files are truncated, but you can click here to view the full file