PageRenderTime 82ms CodeModel.GetById 28ms RepoModel.GetById 1ms app.codeStats 0ms

/wrfv2_fire/phys/module_ra_gsfcsw.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 3138 lines | 1778 code | 504 blank | 856 comment | 46 complexity | 1aade0c232781669363b83e248af8e59 MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. !Comment the following out to turn off aerosol-radiation
    2. 

  2. !feedback between MOSAIC and GSFCSW. wig, 21-Feb-2005
    3. 

  3. 

  4. MODULE module_ra_gsfcsw
    5. 

  5. 

  6.    REAL,    PARAMETER, PRIVATE ::   thresh=1.e-9
    7. 

  7.    REAL,    SAVE               ::   center_lat
    8. 

  8. 

  9. !  Assign co2 and trace gases amount (units are parts/part by volumn)
    10. 

  10. 

  11.    REAL,    PARAMETER, PRIVATE ::   co2   = 300.e-6
    12. 

  12. 

  13. CONTAINS
    14. 

  14. 

  15.    SUBROUTINE GSFCSWRAD(rthraten,gsw,xlat,xlong                   &
    16. 

  16.                    ,dz8w,rho_phy                                  &
    17. 

  17.                    ,alb,t3d,qv3d,qc3d,qr3d                        &
    18. 

  18.                    ,qi3d,qs3d,qg3d,qndrop3d                       &
    19. 

  19.                    ,p3d,p8w3d,pi3d,cldfra3d,rswtoa                &
    20. 

  20.                    ,gmt,cp,g,julday,xtime,declin,solcon           &
    21. 

  21.                    ,radfrq,degrad,taucldi,taucldc,warm_rain       &
    22. 

  22.                    ,tauaer300,tauaer400,tauaer600,tauaer999       & ! jcb
    23. 

  23.                    ,gaer300,gaer400,gaer600,gaer999               & ! jcb
    24. 

  24.                    ,waer300,waer400,waer600,waer999               & ! jcb
    25. 

  25.                    ,aer_ra_feedback                               &
    26. 

  26.                    ,f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop        &
    27. 

  27.                    ,ids,ide, jds,jde, kds,kde                     & 
    28. 

  28.                    ,ims,ime, jms,jme, kms,kme                     &
    29. 

  29.                    ,its,ite, jts,jte, kts,kte                     ) 
    30. 

  30. !------------------------------------------------------------------
    31. 

  31.    IMPLICIT NONE
    32. 

  32. !------------------------------------------------------------------
    33. 

  33.    INTEGER,    PARAMETER     ::        np    = 75
    34. 

  34. 

  35.    INTEGER,    INTENT(IN   ) ::        ids,ide, jds,jde, kds,kde, &
    36. 

  36.                                        ims,ime, jms,jme, kms,kme, &
    37. 

  37.                                        its,ite, jts,jte, kts,kte
    38. 

  38.    LOGICAL,    INTENT(IN   ) ::        warm_rain
    39. 

  39. 

  40.    INTEGER,    INTENT(IN  )  ::                           JULDAY  
    41. 

  41. 

  42. 

  43.    REAL, INTENT(IN    )      ::        RADFRQ,DEGRAD,             &
    44. 

  44.                                        XTIME,DECLIN,SOLCON
    45. 

  45. !
    46. 

  46.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
    47. 

  47.          INTENT(IN    ) ::                                   P3D, &
    48. 

  48.                                                            P8W3D, &
    49. 

  49.                                                             pi3D, &
    50. 

  50.                                                              T3D, &
    51. 

  51.                                                             dz8w, &
    52. 

  52.                                                          rho_phy, &
    53. 

  53.                                                         CLDFRA3D
    54. 

  54. 

  55. 

  56.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
    57. 

  57.          INTENT(INOUT)  ::                              RTHRATEN
    58. 

  58.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
    59. 

  59.          OPTIONAL,                                                &
    60. 

  60.          INTENT(INOUT)  ::                               taucldi, &
    61. 

  61.                                                          taucldc
    62. 

  62. !
    63. 

  63.    REAL, DIMENSION( ims:ime, jms:jme ),                           &
    64. 

  64.          INTENT(IN   )  ::                                  XLAT, &
    65. 

  65.                                                            XLONG, &
    66. 

  66.                                                              ALB
    67. 

  67. !
    68. 

  68.    REAL, DIMENSION( ims:ime, jms:jme ),                           &
    69. 

  69.          INTENT(INOUT)  ::                                   GSW, &
    70. 

  70.                                                           RSWTOA
    71. 

  71. !
    72. 

  72.    REAL, INTENT(IN   )  ::                              GMT,CP,G
    73. 

  73. !
    74. 

  74. 

  75. !
    76. 

  76. ! Optional
    77. 

  77. !
    78. 

  78.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL ,       &
    79. 

  79.          INTENT(IN    ) :: tauaer300,tauaer400,tauaer600,tauaer999, & ! jcb
    80. 

  80.                                  gaer300,gaer400,gaer600,gaer999, & ! jcb
    81. 

  81.                                  waer300,waer400,waer600,waer999    ! jcb
    82. 

  82. 

  83.    INTEGER,    INTENT(IN  ), OPTIONAL   ::       aer_ra_feedback
    84. 

  84. 

  85.    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
    86. 

  86.          OPTIONAL,                                                &           
    87. 

  87.          INTENT(IN    ) ::                                        &
    88. 

  88.                                                             QV3D, &
    89. 

  89.                                                             QC3D, &
    90. 

  90.                                                             QR3D, &
    91. 

  91.                                                             QI3D, &
    92. 

  92.                                                             QS3D, &
    93. 

  93.                                                             QG3D, &
    94. 

  94.                                                         QNDROP3D
    95. 

  95. 

  96.    LOGICAL, OPTIONAL, INTENT(IN ) ::    &
    97. 

  97.                                    F_QV,F_QC,F_QR,F_QI,F_QS,F_QG, &
    98. 

  98.                                                         F_QNDROP
    99. 

  99. 

  100. ! LOCAL VARS
    101. 

  101.  
    102. 

  102.    REAL, DIMENSION( its:ite ) ::                                  &
    103. 

  103.                                                               ts, &
    104. 

  104.                                                             cosz, &
    105. 

  105.                                                               fp, &
    106. 

  106.                                                           rsuvbm, &
    107. 

  107.                                                           rsuvdf, &
    108. 

  108.                                                           rsirbm, &
    109. 

  109.                                                           rsirdf, &
    110. 

  110.                                                             p400, &
    111. 

  111.                                                             p700
    112. 

  112. 

  113.    INTEGER, DIMENSION( its:ite ) ::                               &
    114. 

  114.                                                              ict, &
    115. 

  115.                                                              icb   
    116. 

  116. 

  117.    REAL, DIMENSION( its:ite, kts-1:kte, 2 ) ::            taucld
    118. 

  118. 

  119.    REAL, DIMENSION( its:ite, kts-1:kte+1 )  ::               flx, &
    120. 

  120.                                                             flxd
    121. 

  121. !
    122. 

  122.    REAL, DIMENSION( its:ite, kts-1:kte ) ::                     O3
    123. 

  123. !
    124. 

  124.    REAL, DIMENSION( its:ite, kts-1:kte, 11 ) ::                   &
    125. 

  125.                                                            taual, &
    126. 

  126.                                                            ssaal, &
    127. 

  127.                                                            asyal
    128. 

  128. 

  129.    REAL, DIMENSION( its:ite, kts-1:kte, 2 ) ::                    &
    130. 

  130.                                                             reff, &    
    131. 

  131.                                                              cwc    
    132. 

  132.    REAL, DIMENSION( its: ite, kts-1:kte+1 ) ::                    &
    133. 

  133.                                                            P8W2D
    134. 

  134.    REAL, DIMENSION( its: ite, kts-1:kte ) ::                      &
    135. 

  135.                                                           TTEN2D, &
    136. 

  136.                                                         qndrop2d, &
    137. 

  137.                                                             SH2D, &
    138. 

  138.                                                              P2D, &
    139. 

  139.                                                              T2D, &
    140. 

  140.                                                           fcld2D
    141. 

  141. 

  142.    REAL, DIMENSION( np, 5 ) ::                              pres, &
    143. 

  143.                                                            ozone
    144. 

  144.    REAL, DIMENSION( np )    ::                                 p
    145. 

  145. 

  146.    LOGICAL :: cldwater,overcast, predicate
    147. 

  147. !
    148. 

  148.    INTEGER :: i,j,K,NK,ib,kk,mix,mkx
    149. 

  149. 

  150. !  iprof = 1  :  mid-latitude summer profile
    151. 

  151. !        = 2  :  mid-latitude winter profile
    152. 

  152. !        = 3  :  sub-arctic   summer profile
    153. 

  153. !        = 4  :  sub-arctic   winter profile
    154. 

  154. !        = 5  :  tropical profile
    155. 

  155. !
    156. 

  156. 

  157.    INTEGER  ::                                             iprof, &
    158. 

  158.                                                        is_summer, &
    159. 

  159.                                                        ie_summer, &
    160. 

  160.                                                           lattmp
    161. 

  161. 

  162. 

  163. !
    164. 

  164.    REAL    :: XLAT0,XLONG0
    165. 

  165.    REAL    :: fac,latrmp
    166. 

  166.    REAL    :: xt24,tloctm,hrang,xxlat
    167. 

  167. 

  168.    real, dimension(11) :: midbands  ! jcb
    169. 

  169.    data midbands/.2,.235,.27,.2875,.3025,.305,.3625,.55,1.92,1.745,6.135/   ! jcb
    170. 

  170.    real :: ang,slope ! jcb
    171. 

  171.    character(len=200) :: msg !wig
    172. 

  172.    real pi, third, relconst, lwpmin, rhoh2o
    173. 

  173. !
    174. 

  174. !--------------------------------------------------------------------------------
    175. 

  175. !   data set 1
    176. 

  176. !     mid-latitude summer (75 levels) :  p(mb)  o3(g/g)
    177. 

  177. !     surface temp = 294.0
    178. 

  178. !
    179. 

  179.       data (pres(i,1),i=1,np)/ &
    180. 

  180.           0.0006244,   0.0008759,   0.0012286,   0.0017234,   0.0024174, &
    181. 

  181.           0.0033909,   0.0047565,   0.0066720,   0.0093589,   0.0131278, &
    182. 

  182.           0.0184145,   0.0258302,   0.0362323,   0.0508234,   0.0712906, &
    183. 

  183.           0.1000000,   0.1402710,   0.1967600,   0.2759970,   0.3871430, &
    184. 

  184.           0.5430,    0.7617,    1.0685,    1.4988,    2.1024,    2.9490, &
    185. 

  185.           4.1366,    5.8025,    8.1392,   11.4170,   16.0147,   22.4640, &
    186. 

  186.          31.5105,   44.2001,   62.0000,   85.7750,  109.5500,  133.3250, &
    187. 

  187.         157.1000,  180.8750,  204.6500,  228.4250,  252.2000,  275.9750, &
    188. 

  188.         299.7500,  323.5250,  347.3000,  371.0750,  394.8500,  418.6250, &
    189. 

  189.         442.4000,  466.1750,  489.9500,  513.7250,  537.5000,  561.2750, &
    190. 

  190.         585.0500,  608.8250,  632.6000,  656.3750,  680.1500,  703.9250, &
    191. 

  191.         727.7000,  751.4750,  775.2500,  799.0250,  822.8000,  846.5750, &
    192. 

  192.         870.3500,  894.1250,  917.9000,  941.6750,  965.4500,  989.2250, &
    193. 

  193.        1013.0000/
    194. 

  194. !
    195. 

  195.       data (ozone(i,1),i=1,np)/ &
    196. 

  196.         0.1793E-06,  0.2228E-06,  0.2665E-06,  0.3104E-06,  0.3545E-06, &
    197. 

  197.         0.3989E-06,  0.4435E-06,  0.4883E-06,  0.5333E-06,  0.5786E-06, &
    198. 

  198.         0.6241E-06,  0.6698E-06,  0.7157E-06,  0.7622E-06,  0.8557E-06, &
    199. 

  199.         0.1150E-05,  0.1462E-05,  0.1793E-05,  0.2143E-05,  0.2512E-05, &
    200. 

  200.         0.2902E-05,  0.3313E-05,  0.4016E-05,  0.5193E-05,  0.6698E-05, &
    201. 

  201.         0.8483E-05,  0.9378E-05,  0.9792E-05,  0.1002E-04,  0.1014E-04, &
    202. 

  202.         0.9312E-05,  0.7834E-05,  0.6448E-05,  0.5159E-05,  0.3390E-05, &
    203. 

  203.         0.1937E-05,  0.1205E-05,  0.8778E-06,  0.6935E-06,  0.5112E-06, &
    204. 

  204.         0.3877E-06,  0.3262E-06,  0.2770E-06,  0.2266E-06,  0.2020E-06, &
    205. 

  205.         0.1845E-06,  0.1679E-06,  0.1519E-06,  0.1415E-06,  0.1317E-06, &
    206. 

  206.         0.1225E-06,  0.1137E-06,  0.1055E-06,  0.1001E-06,  0.9487E-07, &
    207. 

  207.         0.9016E-07,  0.8641E-07,  0.8276E-07,  0.7930E-07,  0.7635E-07, &
    208. 

  208.         0.7347E-07,  0.7065E-07,  0.6821E-07,  0.6593E-07,  0.6368E-07, &
    209. 

  209.         0.6148E-07,  0.5998E-07,  0.5859E-07,  0.5720E-07,  0.5582E-07, &
    210. 

  210.         0.5457E-07,  0.5339E-07,  0.5224E-07,  0.5110E-07,  0.4999E-07/
    211. 

  211. 

  212. !--------------------------------------------------------------------------------
    213. 

  213. !   data set 2
    214. 

  214. !   mid-latitude winter (75 levels) :  p(mb)  o3(g/g)
    215. 

  215. !   surface temp = 272.2
    216. 

  216. !
    217. 

  217.       data (pres(i,2),i=1,np)/ &
    218. 

  218.           0.0006244,   0.0008759,   0.0012286,   0.0017234,   0.0024174, &
    219. 

  219.           0.0033909,   0.0047565,   0.0066720,   0.0093589,   0.0131278, &
    220. 

  220.           0.0184145,   0.0258302,   0.0362323,   0.0508234,   0.0712906, &
    221. 

  221.           0.1000000,   0.1402710,   0.1967600,   0.2759970,   0.3871430, &
    222. 

  222.           0.5430,    0.7617,    1.0685,    1.4988,    2.1024,    2.9490, &
    223. 

  223.           4.1366,    5.8025,    8.1392,   11.4170,   16.0147,   22.4640, &
    224. 

  224.          31.5105,   44.2001,   62.0000,   85.9000,  109.8000,  133.7000, &
    225. 

  225.         157.6000,  181.5000,  205.4000,  229.3000,  253.2000,  277.1000, &
    226. 

  226.         301.0000,  324.9000,  348.8000,  372.7000,  396.6000,  420.5000, &
    227. 

  227.         444.4000,  468.3000,  492.2000,  516.1000,  540.0000,  563.9000, &
    228. 

  228.         587.8000,  611.7000,  635.6000,  659.5000,  683.4000,  707.3000, &
    229. 

  229.         731.2000,  755.1000,  779.0000,  802.9000,  826.8000,  850.7000, &
    230. 

  230.         874.6000,  898.5000,  922.4000,  946.3000,  970.2000,  994.1000, &
    231. 

  231.        1018.0000/
    232. 

  232. !
    233. 

  233.       data (ozone(i,2),i=1,np)/ &
    234. 

  234.         0.2353E-06,  0.3054E-06,  0.3771E-06,  0.4498E-06,  0.5236E-06, &
    235. 

  235.         0.5984E-06,  0.6742E-06,  0.7511E-06,  0.8290E-06,  0.9080E-06, &
    236. 

  236.         0.9881E-06,  0.1069E-05,  0.1152E-05,  0.1319E-05,  0.1725E-05, &
    237. 

  237.         0.2145E-05,  0.2581E-05,  0.3031E-05,  0.3497E-05,  0.3980E-05, &
    238. 

  238.         0.4478E-05,  0.5300E-05,  0.6725E-05,  0.8415E-05,  0.1035E-04, &
    239. 

  239.         0.1141E-04,  0.1155E-04,  0.1143E-04,  0.1093E-04,  0.1060E-04, &
    240. 

  240.         0.9720E-05,  0.8849E-05,  0.7424E-05,  0.6023E-05,  0.4310E-05, &
    241. 

  241.         0.2820E-05,  0.1990E-05,  0.1518E-05,  0.1206E-05,  0.9370E-06, &
    242. 

  242.         0.7177E-06,  0.5450E-06,  0.4131E-06,  0.3277E-06,  0.2563E-06, &
    243. 

  243.         0.2120E-06,  0.1711E-06,  0.1524E-06,  0.1344E-06,  0.1199E-06, &
    244. 

  244.         0.1066E-06,  0.9516E-07,  0.8858E-07,  0.8219E-07,  0.7598E-07, &
    245. 

  245.         0.6992E-07,  0.6403E-07,  0.5887E-07,  0.5712E-07,  0.5540E-07, &
    246. 

  246.         0.5370E-07,  0.5214E-07,  0.5069E-07,  0.4926E-07,  0.4785E-07, &
    247. 

  247.         0.4713E-07,  0.4694E-07,  0.4676E-07,  0.4658E-07,  0.4641E-07, &
    248. 

  248.         0.4634E-07,  0.4627E-07,  0.4619E-07,  0.4612E-07,  0.4605E-07/
    249. 

  249. 

  250. 

  251. !--------------------------------------------------------------------------------
    252. 

  252. !   data set 3
    253. 

  253. !   sub-arctic summer (75 levels) :  p(mb)  o3(g/g)
    254. 

  254. !   surface temp = 287.0
    255. 

  255. !
    256. 

  256.       data (pres(i,3),i=1,np)/ &
    257. 

  257.           0.0006244,   0.0008759,   0.0012286,   0.0017234,   0.0024174, &
    258. 

  258.           0.0033909,   0.0047565,   0.0066720,   0.0093589,   0.0131278, &
    259. 

  259.           0.0184145,   0.0258302,   0.0362323,   0.0508234,   0.0712906, &
    260. 

  260.           0.1000000,   0.1402710,   0.1967600,   0.2759970,   0.3871430, &
    261. 

  261.           0.5430,    0.7617,    1.0685,    1.4988,    2.1024,    2.9490, &
    262. 

  262.           4.1366,    5.8025,    8.1392,   11.4170,   16.0147,   22.4640, &
    263. 

  263.          31.5105,   44.2001,   62.0000,   85.7000,  109.4000,  133.1000, &
    264. 

  264.         156.8000,  180.5000,  204.2000,  227.9000,  251.6000,  275.3000, &
    265. 

  265.         299.0000,  322.7000,  346.4000,  370.1000,  393.8000,  417.5000, &
    266. 

  266.         441.2000,  464.9000,  488.6000,  512.3000,  536.0000,  559.7000, &
    267. 

  267.         583.4000,  607.1000,  630.8000,  654.5000,  678.2000,  701.9000, &
    268. 

  268.         725.6000,  749.3000,  773.0000,  796.7000,  820.4000,  844.1000, &
    269. 

  269.         867.8000,  891.5000,  915.2000,  938.9000,  962.6000,  986.3000, &
    270. 

  270.        1010.0000/
    271. 

  271. !
    272. 

  272.       data (ozone(i,3),i=1,np)/ &
    273. 

  273.         0.1728E-06,  0.2131E-06,  0.2537E-06,  0.2944E-06,  0.3353E-06, &
    274. 

  274.         0.3764E-06,  0.4176E-06,  0.4590E-06,  0.5006E-06,  0.5423E-06, &
    275. 

  275.         0.5842E-06,  0.6263E-06,  0.6685E-06,  0.7112E-06,  0.7631E-06, &
    276. 

  276.         0.1040E-05,  0.1340E-05,  0.1660E-05,  0.2001E-05,  0.2362E-05, &
    277. 

  277.         0.2746E-05,  0.3153E-05,  0.3762E-05,  0.4988E-05,  0.6518E-05, &
    278. 

  278.         0.8352E-05,  0.9328E-05,  0.9731E-05,  0.8985E-05,  0.7632E-05, &
    279. 

  279.         0.6814E-05,  0.6384E-05,  0.5718E-05,  0.4728E-05,  0.4136E-05, &
    280. 

  280.         0.3033E-05,  0.2000E-05,  0.1486E-05,  0.1121E-05,  0.8680E-06, &
    281. 

  281.         0.6474E-06,  0.5164E-06,  0.3921E-06,  0.2996E-06,  0.2562E-06, &
    282. 

  282.         0.2139E-06,  0.1723E-06,  0.1460E-06,  0.1360E-06,  0.1267E-06, &
    283. 

  283.         0.1189E-06,  0.1114E-06,  0.1040E-06,  0.9678E-07,  0.8969E-07, &
    284. 

  284.         0.8468E-07,  0.8025E-07,  0.7590E-07,  0.7250E-07,  0.6969E-07, &
    285. 

  285.         0.6694E-07,  0.6429E-07,  0.6208E-07,  0.5991E-07,  0.5778E-07, &
    286. 

  286.         0.5575E-07,  0.5403E-07,  0.5233E-07,  0.5067E-07,  0.4904E-07, &
    287. 

  287.         0.4721E-07,  0.4535E-07,  0.4353E-07,  0.4173E-07,  0.3997E-07/
    288. 

  288. 

  289. 

  290. !--------------------------------------------------------------------------------
    291. 

  291. !   data set 3
    292. 

  292. !   sub-arctic winter (75 levels) :   p(mb)  o3(g/g)
    293. 

  293. !   surface temp = 257.1
    294. 

  294. !
    295. 

  295.       data (pres(i,4),i=1,np)/ &
    296. 

  296.           0.0006244,   0.0008759,   0.0012286,   0.0017234,   0.0024174, &
    297. 

  297.           0.0033909,   0.0047565,   0.0066720,   0.0093589,   0.0131278, &
    298. 

  298.           0.0184145,   0.0258302,   0.0362323,   0.0508234,   0.0712906, &
    299. 

  299.           0.1000000,   0.1402710,   0.1967600,   0.2759970,   0.3871430, &
    300. 

  300.           0.5430,    0.7617,    1.0685,    1.4988,    2.1024,    2.9490, &
    301. 

  301.           4.1366,    5.8025,    8.1392,   11.4170,   16.0147,   22.4640, &
    302. 

  302.          31.5105,   44.2001,   62.0000,   85.7750,  109.5500,  133.3250, &
    303. 

  303.         157.1000,  180.8750,  204.6500,  228.4250,  252.2000,  275.9750, &
    304. 

  304.         299.7500,  323.5250,  347.3000,  371.0750,  394.8500,  418.6250, &
    305. 

  305.         442.4000,  466.1750,  489.9500,  513.7250,  537.5000,  561.2750, &
    306. 

  306.         585.0500,  608.8250,  632.6000,  656.3750,  680.1500,  703.9250, &
    307. 

  307.         727.7000,  751.4750,  775.2500,  799.0250,  822.8000,  846.5750, &
    308. 

  308.         870.3500,  894.1250,  917.9000,  941.6750,  965.4500,  989.2250, &
    309. 

  309.        1013.0000/
    310. 

  310. !
    311. 

  311.       data (ozone(i,4),i=1,np)/ &
    312. 

  312.         0.2683E-06,  0.3562E-06,  0.4464E-06,  0.5387E-06,  0.6333E-06, &
    313. 

  313.         0.7301E-06,  0.8291E-06,  0.9306E-06,  0.1034E-05,  0.1140E-05, &
    314. 

  314.         0.1249E-05,  0.1360E-05,  0.1474E-05,  0.1855E-05,  0.2357E-05, &
    315. 

  315.         0.2866E-05,  0.3383E-05,  0.3906E-05,  0.4437E-05,  0.4975E-05, &
    316. 

  316.         0.5513E-05,  0.6815E-05,  0.8157E-05,  0.1008E-04,  0.1200E-04, &
    317. 

  317.         0.1242E-04,  0.1250E-04,  0.1157E-04,  0.1010E-04,  0.9063E-05, &
    318. 

  318.         0.8836E-05,  0.8632E-05,  0.8391E-05,  0.7224E-05,  0.6054E-05, &
    319. 

  319.         0.4503E-05,  0.3204E-05,  0.2278E-05,  0.1833E-05,  0.1433E-05, &
    320. 

  320.         0.9996E-06,  0.7440E-06,  0.5471E-06,  0.3944E-06,  0.2852E-06, &
    321. 

  321.         0.1977E-06,  0.1559E-06,  0.1333E-06,  0.1126E-06,  0.9441E-07, &
    322. 

  322.         0.7678E-07,  0.7054E-07,  0.6684E-07,  0.6323E-07,  0.6028E-07, &
    323. 

  323.         0.5746E-07,  0.5468E-07,  0.5227E-07,  0.5006E-07,  0.4789E-07, &
    324. 

  324.         0.4576E-07,  0.4402E-07,  0.4230E-07,  0.4062E-07,  0.3897E-07, &
    325. 

  325.         0.3793E-07,  0.3697E-07,  0.3602E-07,  0.3506E-07,  0.3413E-07, &
    326. 

  326.         0.3326E-07,  0.3239E-07,  0.3153E-07,  0.3069E-07,  0.2987E-07/ 
    327. 

  327. 

  328. !--------------------------------------------------------------------------------
    329. 

  329. !   data set 4
    330. 

  330. !   tropical (75 levels) :   p(mb)  o3(g/g)
    331. 

  331. !   surface temp = 300.0
    332. 

  332. !
    333. 

  333.       data (pres(i,5),i=1,np)/ &
    334. 

  334.           0.0006244,   0.0008759,   0.0012286,   0.0017234,   0.0024174, &
    335. 

  335.           0.0033909,   0.0047565,   0.0066720,   0.0093589,   0.0131278, &
    336. 

  336.           0.0184145,   0.0258302,   0.0362323,   0.0508234,   0.0712906, &
    337. 

  337.           0.1000000,   0.1402710,   0.1967600,   0.2759970,   0.3871430, &
    338. 

  338.           0.5430,    0.7617,    1.0685,    1.4988,    2.1024,    2.9490, &
    339. 

  339.           4.1366,    5.8025,    8.1392,   11.4170,   16.0147,   22.4640, &
    340. 

  340.          31.5105,   44.2001,   62.0000,   85.7750,  109.5500,  133.3250, &
    341. 

  341.         157.1000,  180.8750,  204.6500,  228.4250,  252.2000,  275.9750, &
    342. 

  342.         299.7500,  323.5250,  347.3000,  371.0750,  394.8500,  418.6250, &
    343. 

  343.         442.4000,  466.1750,  489.9500,  513.7250,  537.5000,  561.2750, &
    344. 

  344.         585.0500,  608.8250,  632.6000,  656.3750,  680.1500,  703.9250, &
    345. 

  345.         727.7000,  751.4750,  775.2500,  799.0250,  822.8000,  846.5750, &
    346. 

  346.         870.3500,  894.1250,  917.9000,  941.6750,  965.4500,  989.2250, &
    347. 

  347.        1013.0000/
    348. 

  348. !
    349. 

  349.       data (ozone(i,5),i=1,np)/ &
    350. 

  350.         0.1993E-06,  0.2521E-06,  0.3051E-06,  0.3585E-06,  0.4121E-06, &
    351. 

  351.         0.4661E-06,  0.5203E-06,  0.5748E-06,  0.6296E-06,  0.6847E-06, &
    352. 

  352.         0.7402E-06,  0.7959E-06,  0.8519E-06,  0.9096E-06,  0.1125E-05, &
    353. 

  353.         0.1450E-05,  0.1794E-05,  0.2156E-05,  0.2538E-05,  0.2939E-05, &
    354. 

  354.         0.3362E-05,  0.3785E-05,  0.4753E-05,  0.6005E-05,  0.7804E-05, &
    355. 

  355.         0.9635E-05,  0.1023E-04,  0.1067E-04,  0.1177E-04,  0.1290E-04, &
    356. 

  356.         0.1134E-04,  0.9223E-05,  0.6667E-05,  0.3644E-05,  0.1545E-05, &
    357. 

  357.         0.5355E-06,  0.2523E-06,  0.2062E-06,  0.1734E-06,  0.1548E-06, &
    358. 

  358.         0.1360E-06,  0.1204E-06,  0.1074E-06,  0.9707E-07,  0.8960E-07, &
    359. 

  359.         0.8419E-07,  0.7962E-07,  0.7542E-07,  0.7290E-07,  0.7109E-07, &
    360. 

  360.         0.6940E-07,  0.6786E-07,  0.6635E-07,  0.6500E-07,  0.6370E-07, &
    361. 

  361.         0.6244E-07,  0.6132E-07,  0.6022E-07,  0.5914E-07,  0.5884E-07, &
    362. 

  362.         0.5855E-07,  0.5823E-07,  0.5772E-07,  0.5703E-07,  0.5635E-07, &
    363. 

  363.         0.5570E-07,  0.5492E-07,  0.5412E-07,  0.5335E-07,  0.5260E-07, &
    364. 

  364.         0.5167E-07,  0.5063E-07,  0.4961E-07,  0.4860E-07,  0.4761E-07/
    365. 

  365. 

  366. !--------------------------------------------------------------------------------
    367. 

  367. 

  368. #ifdef WRF_CHEM
    369. 

  369.    IF ( aer_ra_feedback == 1) then
    370. 

  370.    IF ( .NOT. &
    371. 

  371.       ( PRESENT(tauaer300) .AND. &
    372. 

  372.         PRESENT(tauaer400) .AND. &
    373. 

  373.         PRESENT(tauaer600) .AND. &
    374. 

  374.         PRESENT(tauaer999) .AND. &
    375. 

  375.         PRESENT(gaer300) .AND. &
    376. 

  376.         PRESENT(gaer400) .AND. &
    377. 

  377.         PRESENT(gaer600) .AND. &
    378. 

  378.         PRESENT(gaer999) .AND. &
    379. 

  379.         PRESENT(waer300) .AND. &
    380. 

  380.         PRESENT(waer400) .AND. &
    381. 

  381.         PRESENT(waer600) .AND. &
    382. 

  382.         PRESENT(waer999) ) ) THEN
    383. 

  383.       CALL wrf_error_fatal ( 'Warning: missing fields required for aerosol radiation' )
    384. 

  384.    ENDIF
    385. 

  385.    ENDIF
    386. 

  386. #endif
    387. 

  387.    cldwater = .true.
    388. 

  388.    overcast = .false.
    389. 

  389. 

  390.    mix=ite-its+1 
    391. 

  391.    mkx=kte-kts+1 
    392. 

  392. 

  393.    is_summer=80
    394. 

  394.    ie_summer=265
    395. 

  395. 

  396. ! testing, need to change iprof, which is function of lat and julian day
    397. 

  397. !  iprof = 1  :  mid-latitude summer profile
    398. 

  398. !        = 2  :  mid-latitude winter profile
    399. 

  399. !        = 3  :  sub-arctic   summer profile
    400. 

  400. !        = 4  :  sub-arctic   winter profile
    401. 

  401. !        = 5  :  tropical profile
    402. 

  402. 

  403.    IF (abs(center_lat) .le. 30. ) THEN ! tropic
    404. 

  404.       iprof = 5
    405. 

  405.    ELSE
    406. 

  406.       IF (center_lat .gt.  0.) THEN
    407. 

  407.          IF (center_lat .gt. 60. ) THEN !  arctic
    408. 

  408.             IF (JULDAY .gt. is_summer .and. JULDAY .lt. ie_summer ) THEN
    409. 

  409.                ! arctic summer
    410. 

  410.                iprof = 3
    411. 

  411.             ELSE
    412. 

  412.                ! arctic winter
    413. 

  413.                iprof = 4
    414. 

  414.             ENDIF
    415. 

  415.          ELSE        ! midlatitude
    416. 

  416.             IF (JULDAY .gt. is_summer .and. JULDAY .lt. ie_summer ) THEN
    417. 

  417.                ! north midlatitude summer
    418. 

  418.                iprof = 1
    419. 

  419.             ELSE
    420. 

  420.                ! north midlatitude winter
    421. 

  421.                iprof = 2
    422. 

  422.             ENDIF
    423. 

  423.          ENDIF
    424. 

  424. 

  425.       ELSE
    426. 

  426.          IF (center_lat .lt. -60. ) THEN !  antarctic
    427. 

  427.             IF (JULDAY .lt. is_summer .or. JULDAY .gt. ie_summer ) THEN
    428. 

  428.                ! antarctic summer
    429. 

  429.                iprof = 3
    430. 

  430.             ELSE
    431. 

  431.                ! antarctic winter
    432. 

  432.                iprof = 4
    433. 

  433.             ENDIF
    434. 

  434.          ELSE        ! midlatitude
    435. 

  435.             IF (JULDAY .lt. is_summer .or. JULDAY .gt. ie_summer ) THEN
    436. 

  436.                ! south midlatitude summer
    437. 

  437.                iprof = 1
    438. 

  438.             ELSE
    439. 

  439.                ! south midlatitude winter
    440. 

  440.                iprof = 2
    441. 

  441.             ENDIF
    442. 

  442.          ENDIF
    443. 

  443. 

  444.       ENDIF
    445. 

  445.    ENDIF
    446. 

  446. 

  447. 

  448.    j_loop: DO J=jts,jte
    449. 

  449. 

  450.       DO K=kts,kte          
    451. 

  451.       DO I=its,ite          
    452. 

  452.          cwc(i,k,1) = 0.
    453. 

  453.          cwc(i,k,2) = 0.
    454. 

  454.       ENDDO
    455. 

  455.       ENDDO
    456. 

  456. 

  457.       DO K=1,np
    458. 

  458.          p(k)=pres(k,iprof)
    459. 

  459.       ENDDO
    460. 

  460. 

  461. ! reverse vars 
    462. 

  462. !
    463. 

  463.       DO K=kts,kte+1
    464. 

  464.       DO I=its,ite
    465. 

  465.          NK=kme-K+kms
    466. 

  466.          P8W2D(I,K)=p8w3d(i,nk,j)*0.01   ! P8w2D is in mb
    467. 

  467.       ENDDO
    468. 

  468.       ENDDO
    469. 

  469. 

  470.       DO I=its,ite
    471. 

  471.          P8W2D(I,0)=.0
    472. 

  472.       ENDDO
    473. 

  473. !
    474. 

  474.       DO K=kts,kte
    475. 

  475.       DO I=its,ite
    476. 

  476.          NK=kme-1-K+kms
    477. 

  477.          TTEN2D(I,K)=0.
    478. 

  478.          T2D(I,K)=T3D(I,NK,J)
    479. 

  479. 

  480. ! SH2D specific humidity
    481. 

  481.          SH2D(I,K)=QV3D(I,NK,J)/(1.+QV3D(I,NK,J))
    482. 

  482.          SH2D(I,K)=max(0.,SH2D(I,K))
    483. 

  483.          cwc(I,K,2)=QC3D(I,NK,J)
    484. 

  484.          cwc(I,K,2)=max(0.,cwc(I,K,2))
    485. 

  485. 

  486.          P2D(I,K)=p3d(i,nk,j)*0.01      ! P2D is in mb
    487. 

  487.          fcld2D(I,K)=CLDFRA3D(I,NK,J)
    488. 

  488.       ENDDO
    489. 

  489.       ENDDO
    490. 

  490. 

  491. ! This logic is tortured because cannot test F_QI unless
    492. 

  492. ! it is present, and order of evaluation of expressions
    493. 

  493. ! is not specified in Fortran
    494. 

  494. 

  495.       IF ( PRESENT ( F_QI ) ) THEN
    496. 

  496.          predicate = F_QI
    497. 

  497.       ELSE
    498. 

  498.         predicate = .FALSE.
    499. 

  499.       ENDIF
    500. 

  500. 

  501.       IF (.NOT. warm_rain .AND. .NOT. predicate ) THEN
    502. 

  502.          DO K=kts,kte
    503. 

  503.          DO I=its,ite
    504. 

  504.             IF (T2D(I,K) .lt. 273.15) THEN
    505. 

  505.                cwc(I,K,1)=cwc(I,K,2)
    506. 

  506.                cwc(I,K,2)=0.
    507. 

  507.             ENDIF
    508. 

  508.          ENDDO
    509. 

  509.          ENDDO
    510. 

  510.       ENDIF
    511. 

  511. 

  512.        IF ( PRESENT( F_QNDROP ) ) THEN
    513. 

  513.          IF ( F_QNDROP ) THEN
    514. 

  514.             DO K=kts,kte
    515. 

  515.                DO I=its,ite
    516. 

  516.                   NK=kme-1-K+kms
    517. 

  517.                   qndrop2d(I,K)=qndrop3d(I,NK,j)
    518. 

  518.                ENDDO
    519. 

  519.             ENDDO
    520. 

  520.             qndrop2d(:,kts-1)=0.
    521. 

  521.          END IF
    522. 

  522.       END IF
    523. 

  523. 

  524.      DO I=its,ite
    525. 

  525.          TTEN2D(I,0)=0.
    526. 

  526.          T2D(I,0)=T2D(I,1)
    527. 

  527. ! SH2D specific humidity
    528. 

  528.          SH2D(I,0)=0.5*SH2D(i,1)
    529. 

  529.          cwc(I,0,2)=0.
    530. 

  530.          cwc(I,0,1)=0.
    531. 

  531.          P2D(I,0)=0.5*(P8W2D(I,0)+P8W2D(I,1))
    532. 

  532.          fcld2D(I,0)=0.
    533. 

  533.       ENDDO
    534. 

  534. !
    535. 

  535.       IF ( PRESENT( F_QI ) .AND. PRESENT( qi3d)  ) THEN
    536. 

  536. 	 IF ( (F_QI) ) THEN
    537. 

  537.             DO K=kts,kte          
    538. 

  538.             DO I=its,ite          
    539. 

  539.                NK=kme-1-K+kms
    540. 

  540.                cwc(I,K,1)=QI3D(I,NK,J)
    541. 

  541.                cwc(I,K,1)=max(0.,cwc(I,K,1))
    542. 

  542.             ENDDO
    543. 

  543.             ENDDO
    544. 

  544.          ENDIF
    545. 

  545.       ENDIF
    546. 

  546. !
    547. 

  547. ! ... Vertical profiles for ozone
    548. 

  548. !
    549. 

  549.       call o3prof (np, p, ozone(1,iprof), its, ite, kts-1, kte, P2D, O3)
    550. 

  550. 

  551. ! ... Vertical profiles for effective particle size
    552. 

  552. !
    553. 

  553.       pi = 4.*atan(1.0)
    554. 

  554.       third=1./3.
    555. 

  555.       rhoh2o=1.e3
    556. 

  556.       relconst=3/(4.*pi*rhoh2o)
    557. 

  557. !     minimun liquid water path to calculate rel
    558. 

  558. !     corresponds to optical depth of 1.e-3 for radius 4 microns.
    559. 

  559.       lwpmin=3.e-5
    560. 

  560.       do k = kts-1, kte
    561. 

  561.       do i = its, ite
    562. 

  562.          reff(i,k,2) = 10.
    563. 

  563.          if( PRESENT( F_QNDROP ) ) then
    564. 

  564.             if( F_QNDROP ) then
    565. 

  565.               if ( cwc(i,k,2)*(P8W2D(I,K+1)-P8W2D(I,K)).gt.lwpmin.and. &
    566. 

  566.                    qndrop2d(i,k).gt.1000. ) then
    567. 

  567.                reff(i,k,2)=(relconst*cwc(i,k,2)/qndrop2d(i,k))**third ! effective radius in m
    568. 

  568. !           apply scaling from Martin et al., JAS 51, 1830.
    569. 

  569.                reff(i,k,2)=1.1*reff(i,k,2)
    570. 

  570.                reff(i,k,2)=reff(i,k,2)*1.e6 ! convert from m to microns
    571. 

  571.                reff(i,k,2)=max(reff(i,k,2),4.)
    572. 

  572.                reff(i,k,2)=min(reff(i,k,2),20.)
    573. 

  573.               end if
    574. 

  574.             end if
    575. 

  575.          end if
    576. 

  576.          reff(i,k,1) = 80.
    577. 

  577.       end do
    578. 

  578.       end do
    579. 

  579. !
    580. 

  580. ! ... Level indices separating high, middle and low clouds
    581. 

  581. !
    582. 

  582.       do i = its, ite
    583. 

  583.          p400(i) = 1.e5
    584. 

  584.          p700(i) = 1.e5
    585. 

  585.       enddo
    586. 

  586. 

  587.       do k = kts-1,kte+1
    588. 

  588.          do i = its, ite
    589. 

  589.             if (abs(P8W2D(i,k) - 400.) .lt. p400(i)) then
    590. 

  590.                p400(i) = abs(P8W2D(i,k) - 400.)
    591. 

  591.                ict(i) = k
    592. 

  592.             endif
    593. 

  593.             if (abs(P8W2D(i,k) - 700.) .lt. p700(i)) then
    594. 

  594.                p700(i) = abs(P8W2D(i,k) - 700.)
    595. 

  595.                icb(i) = k
    596. 

  596.             endif
    597. 

  597.         end do
    598. 

  598.       end do
    599. 

  599. 

  600. !wig beg
    601. 

  601. ! ... Aerosol effects. Added aerosol feedbacks with MOSAIC, Dec. 2005.
    602. 

  602. !
    603. 

  603.       do ib = 1, 11
    604. 

  604.       do k = kts-1,kte
    605. 

  605.       do i = its,ite
    606. 

  606.          taual(i,k,ib) = 0.
    607. 

  607.          ssaal(i,k,ib) = 0.
    608. 

  608.          asyal(i,k,ib) = 0.
    609. 

  609.       end do
    610. 

  610.       end do
    611. 

  611.       end do
    612. 

  612. 

  613. #ifdef WRF_CHEM
    614. 

  614.    IF ( AER_RA_FEEDBACK == 1) then
    615. 

  615. !wig end
    616. 

  616.       do ib = 1, 11
    617. 

  617.       do k = kts-1,kte-1      !wig
    618. 

  618.       do i = its,ite
    619. 

  619. 

  620. !        taual(i,kte-k,ib) = 0.
    621. 

  621. !        ssaal(i,kte-k,ib) = 0.
    622. 

  622. !        asyal(i,kte-k,ib) = 0.
    623. 

  623. 

  624. !jcb beg
    625. 

  625. ! convert optical properties at 300,400,600, and 999 to conform to the band wavelengths
    626. 

  626. ! these are: 200,235,270,287.5,302.5,305,362.5,550,1920,1745,6135; why the emphasis on the UV?
    627. 

  627. ! taual - use angstrom exponent
    628. 

  628.         if(tauaer300(i,k+1,j).gt.thresh .and. tauaer999(i,k+1,j).gt.thresh) then
    629. 

  629.            ang=log(tauaer300(i,k+1,j)/tauaer999(i,k+1,j))/log(999./300.)
    630. 

  630. !       write(6,*)i,k,ang,tauaer300(i,k+1,j),tauaer999(i,k+1,j)
    631. 

  631.            taual(i,kte-k,ib)=tauaer400(i,k+1,j)*(0.4/midbands(ib))**ang ! notice reserved variable
    632. 

  632. !       write(6,10001)i,k,ang,tauaer300(i,k+1,j),tauaer999(i,k+1,j),midbands(ib),taual(i,k,ib)
    633. 

  633. !10001      format(i3,i3,5f12.6)
    634. 

  634. 

  635. ! ssa - linear interpolation; extrapolation
    636. 

  636.            slope=(waer600(i,k+1,j)-waer400(i,k+1,j))/.2
    637. 

  637.            ssaal(i,kte-k,ib) = slope*(midbands(ib)-.6)+waer600(i,k+1,j) ! notice reversed variables
    638. 

  638.            if(ssaal(i,kte-k,ib).lt.0.4) ssaal(i,kte-k,ib)=0.4
    639. 

  639.            if(ssaal(i,kte-k,ib).ge.1.0) ssaal(i,kte-k,ib)=1.0
    640. 

  640. 

  641. ! g - linear interpolation;extrapolation
    642. 

  642.            slope=(gaer600(i,k+1,j)-gaer400(i,k+1,j))/.2
    643. 

  643.            asyal(i,kte-k,ib) = slope*(midbands(ib)-.6)+gaer600(i,k+1,j) ! notice reversed varaibles
    644. 

  644.            if(asyal(i,kte-k,ib).lt.0.5) asyal(i,kte-k,ib)=0.5
    645. 

  645.            if(asyal(i,kte-k,ib).ge.1.0) asyal(i,kte-k,ib)=1.0
    646. 

  646.         endif
    647. 

  647. !jcb end
    648. 

  648.       end do
    649. 

  649.       end do
    650. 

  650.       end do
    651. 

  651. 

  652. !wig beg
    653. 

  653.       do ib = 1, 11
    654. 

  654.       do i = its,ite
    655. 

  655.          slope = 0.  !use slope as a sum holder
    656. 

  656.          do k = kts-1,kte
    657. 

  657.             slope = slope + taual(i,k,ib)
    658. 

  658.          end do
    659. 

  659.          if( slope < 0. ) then
    660. 

  660.             write(msg,'("ERROR: Negative total optical depth of ",f8.2," at point i,j,ib=",3i5)') slope,i,j,ib
    661. 

  661.             call wrf_error_fatal(msg)
    662. 

  662.          else if( slope > 5. ) then
    663. 

  663.             call wrf_message("-------------------------")
    664. 

  664.             write(msg,'("WARNING: Large total optical depth of ",f8.2," at point i,j,ib=",3i5)') slope,i,j,ib
    665. 

  665.             call wrf_message(msg)
    666. 

  666. 

  667.             call wrf_message("Diagnostics 1: k, tauaer300, tauaer400, tauaer600, tauaer999")
    668. 

  668.             do k=kts,kte
    669. 

  669.                write(msg,'(i4,4f8.2)') k, tauaer300(i,k,j), tauaer400(i,k,j), &
    670. 

  670.                     tauaer600(i,k,j), tauaer999(i,k,j)
    671. 

  671.                call wrf_message(msg)
    672. 

  672.             end do
    673. 

  673. 

  674.             call wrf_message("Diagnostics 2: k, gaer300, gaer400, gaer600, gaer999")
    675. 

  675.             do k=kts,kte
    676. 

  676.                write(msg,'(i4,4f8.2)') k, gaer300(i,k,j), gaer400(i,k,j), &
    677. 

  677.                     gaer600(i,k,j), gaer999(i,k,j)
    678. 

  678.                call wrf_message(msg)
    679. 

  679.             end do
    680. 

  680. 

  681.             call wrf_message("Diagnostics 3: k, waer300, waer400, waer600, waer999")
    682. 

  682.             do k=kts,kte
    683. 

  683.                write(msg,'(i4,4f8.2)') k, waer300(i,k,j), waer400(i,k,j), &
    684. 

  684.                     waer600(i,k,j), waer999(i,k,j)
    685. 

  685.                call wrf_message(msg)
    686. 

  686.             end do
    687. 

  687. 

  688.             call wrf_message("Diagnostics 4: k, ssaal, asyal, taual")
    689. 

  689.             do k=kts-1,kte
    690. 

  690.                write(msg,'(i4,3f8.2)') k, ssaal(i,k,ib), asyal(i,k,ib), taual(i,k,ib)
    691. 

  691.                call wrf_message(msg)
    692. 

  692.             end do
    693. 

  693.             call wrf_message("-------------------------")
    694. 

  694.          end if
    695. 

  695.       end do
    696. 

  696.       end do
    697. 

  697. !wig end
    698. 

  698.       endif
    699. 

  699. #endif
    700. 

  700. !
    701. 

  701. ! ... Initialize output arrays
    702. 

  702. !
    703. 

  703.       do ib = 1, 2
    704. 

  704.       do k = kts-1, kte
    705. 

  705.       do i = its, ite
    706. 

  706.          taucld(i,k,ib) = 0.
    707. 

  707.       end do
    708. 

  708.       end do
    709. 

  709.       end do
    710. 

  710. !
    711. 

  711.       do k = kts-1,kte+1
    712. 

  712.       do i = its,ite
    713. 

  713.          flx(i,k)   = 0.
    714. 

  714.          flxd(i,k)  = 0.
    715. 

  715.       end do
    716. 

  716.       end do
    717. 

  717. !
    718. 

  718. ! ... Solar zenith angle
    719. 

  719. !
    720. 

  720.       do i = its,ite
    721. 

  721.         xt24 = mod(xtime + radfrq * 0.5, 1440.)
    722. 

  722.         tloctm = GMT + xt24 / 60. + XLONG(i,j) / 15.
    723. 

  723.         hrang = 15. * (tloctm - 12.) * degrad
    724. 

  724.         xxlat = XLAT(i,j) * degrad
    725. 

  725.         cosz(i) = sin(xxlat) * sin(declin) + &
    726. 

  726.                   cos(xxlat) * cos(declin) * cos(hrang)
    727. 

  727.         rsuvbm(i) = ALB(i,j)
    728. 

  728.         rsuvdf(i) = ALB(i,j)
    729. 

  729.         rsirbm(i) = ALB(i,j)
    730. 

  730.         rsirdf(i) = ALB(i,j)
    731. 

  731.       end do
    732. 

  732.                                   
    733. 

  733.       call sorad (mix,1,1,mkx+1,p8w2D,t2D,sh2D,o3,                 &
    734. 

  734.                   overcast,cldwater,cwc,taucld,reff,fcld2D,ict,icb,&
    735. 

  735.                   taual,ssaal,asyal,                               &
    736. 

  736.                   cosz,rsuvbm,rsuvdf,rsirbm,rsirdf,                &
    737. 

  737.                   flx,flxd)
    738. 

  738. !
    739. 

  739. ! ... Convert the units of flx and flc from fraction to w/m^2
    740. 

  740. !
    741. 

  741.       do k = kts, kte
    742. 

  742.       do i = its, ite
    743. 

  743.          nk=kme-1-k+kms
    744. 

  744.          if(present(taucldc)) taucldc(i,nk,j)=taucld(i,k,2)
    745. 

  745.          if(present(taucldi)) taucldi(i,nk,j)=taucld(i,k,1)
    746. 

  746.       enddo
    747. 

  747.       enddo
    748. 

  748.  
    749. 

  749.       do k = kts, kte+1
    750. 

  750.         do i = its, ite
    751. 

  751.           if (cosz(i) .lt. thresh) then
    752. 

  752.             flx(i,k) = 0.
    753. 

  753.           else
    754. 

  754.             flx(i,k) = flx(i,k) * SOLCON * cosz(i)
    755. 

  755.           endif
    756. 

  756.         end do
    757. 

  757.       end do
    758. 

  758. !
    759. 

  759. ! ... Calculate heating rate (deg/sec)
    760. 

  760. !
    761. 

  761.       fac = .01 * g / Cp
    762. 

  762.       do k = kts, kte
    763. 

  763.       do i = its, ite
    764. 

  764.          if (cosz(i) .gt. thresh) then
    765. 

  765.              TTEN2D(i,k) = - fac * (flx(i,k) - flx(i,k+1))/ &
    766. 

  766.                            (p8w2d(i,k)-p8w2d(i,k+1))
    767. 

  767.          endif
    768. 

  768.       end do
    769. 

  769.       end do
    770. 

  770. 

  771. !     upward top of atmosphere
    772. 

  772.       do i = its, ite
    773. 

  773.         if (cosz(i) .le. thresh) then
    774. 

  774.           RSWTOA(i,j) = 0.
    775. 

  775.         else
    776. 

  776.           RSWTOA(i,j) = flx(i,kts) - flxd(i,kts) * SOLCON * cosz(i)
    777. 

  777.         endif
    778. 

  778.       end do
    779. 

  779. !
    780. 

  780. ! ... Absorbed part in surface energy budget
    781. 

  781. !
    782. 

  782.       do i = its, ite
    783. 

  783.         if (cosz(i) .le. thresh) then
    784. 

  784.           GSW(i,j) = 0.
    785. 

  785.         else
    786. 

  786.           GSW(i,j) = (1. - rsuvbm(i)) * flxd(i,kte+1) * SOLCON * cosz(i)
    787. 

  787.         endif
    788. 

  788.       end do
    789. 

  789. 

  790.       DO K=kts,kte          
    791. 

  791.          NK=kme-1-K+kms
    792. 

  792.          DO I=its,ite
    793. 

  793. ! FIX FROM GODDARD FOR NEGATIVE VALUES
    794. 

  794.             TTEN2D(I,NK)=MAX(TTEN2D(I,NK),0.)
    795. 

  795.             RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+TTEN2D(I,NK)/pi3D(I,K,J)
    796. 

  796.          ENDDO
    797. 

  797.       ENDDO
    798. 

  798. !
    799. 

  799.    ENDDO j_loop                                          
    800. 

  800. 

  801.    END SUBROUTINE GSFCSWRAD
    802. 

  802. 

  803. !*********************   Version Solar-6 (May 8, 1997)  *****************
    804. 

  804. 

  805.       subroutine sorad (m,n,ndim,np,pl,ta,wa,oa,                        &
    806. 

  806.                         overcast,cldwater,cwc,taucld,reff,fcld,ict,icb, &
    807. 

  807.                         taual,ssaal,asyal,                              &
    808. 

  808.                         cosz,rsuvbm,rsuvdf,rsirbm,rsirdf,               &
    809. 

  809.                         flx,flxd)
    810. 

  810. 

  811. !************************************************************************
    812. 

  812. !
    813. 

  813. !                        Version Solar-6 (May 8, 1997)
    814. 

  814. !
    815. 

  815. !  New feature of this version is:
    816. 

  816. !   (1) An option is added for scaling the cloud optical thickness. If
    817. 

  817. !       the fractional cloud cover, fcld, in an atmospheric model is alway 
    818. 

  818. !       either 1 or 0 (i.e. partly cloudy sky is not allowed), it does
    819. 

  819. !       not require the scaling of cloud optical thickness, and the
    820. 

  820. !       option "overcast" can be set to .true.  Computation is faster
    821. 

  821. !       with this option than with overcast=.false.
    822. 

  822. !
    823. 

  823. !**********************************************************************
    824. 

  824. !
    825. 

  825. !                        Version Solar-5 (April 1997)
    826. 

  826. !
    827. 

  827. !  New features of this version are:
    828. 

  828. !   (1) Cloud optical properties can be computed from cloud water/ice
    829. 

  829. !       amount and the effective particle size.
    830. 

  830. !   (2) Aerosol optical properties are functions of height and band.
    831. 

  831. !   (3) A maximum-random cloud overlapping approximation is applied.
    832. 

  832. !
    833. 

  833. !*********************************************************************
    834. 

  834. !  
    835. 

  835. ! This routine computes solar fluxes due to the absoption by water
    836. 

  836. !  vapor, ozone, co2, o2, clouds, and aerosols and due to the
    837. 

  837. !  scattering by clouds, aerosols, and gases.
    838. 

  838. !
    839. 

  839. ! The solar spectrum is divided into one UV+visible band and three IR
    840. 

  840. !  bands separated by the wavelength 0.7 micron.  The UV+visible band
    841. 

  841. !  is further divided into eight sub-bands.
    842. 

  842. !
    843. 

  843. ! This is a vectorized code. It computes fluxes simultaneously for
    844. 

  844. !  (m x n) soundings, which is a subset of (m x ndim) soundings.
    845. 

  845. !  In a global climate model, m and ndim correspond to the numbers of
    846. 

  846. !  grid boxes in the zonal and meridional directions, respectively.
    847. 

  847. !
    848. 

  848. ! Ice and liquid cloud particles are allowed to co-exist in a layer. 
    849. 

  849. !
    850. 

  850. ! There is an option of providing either cloud ice/water mixing ratio 
    851. 

  851. !  (cwc) or thickness (taucld).  If the former is provided, set
    852. 

  852. !  cldwater=.true., and taucld will be computed from cwc and reff as a
    853. 

  853. !  function of spectra band. Otherwise, set cldwater=.false., and
    854. 

  854. !  specify taucld, independent of spectral band.
    855. 

  855. !
    856. 

  856. ! If no information is available for reff, a default value of
    857. 

  857. !  10 micron for liquid water and 75 micron for ice can be used.
    858. 

  858. !  For a clear layer, reff can be set to any values except zero.
    859. 

  859. !
    860. 

  860. ! The maximum-random assumption is applied for treating cloud
    861. 

  861. !  overlapping.
    862. 

  862. 

  863. ! Clouds are grouped into high, middle, and low clouds separated by
    864. 

  864. !  the level indices ict and icb.  For detail, see subroutine cldscale.
    865. 

  865. !
    866. 

  866. ! In a high spatial-resolution atmospheric model, fractional cloud cover
    867. 

  867. !  might be computed to be either 0 or 1.  In such a case, scaling of the
    868. 

  868. !  cloud optical thickness is not necessary, and the computation can be
    869. 

  869. !  made faster by setting overcast=.true.  The option overcast=.false.
    870. 

  870. !  can be applied to any values of the fractional cloud cover, but the
    871. 

  871. !  computation is slower.
    872. 

  872. !
    873. 

  873. ! Aerosol optical thickness, single-scattering albaedo, and asymmtry
    874. 

  874. !  factor can be specified as functions of height and spectral band.
    875. 

  875. !
    876. 

  876. !----- Input parameters:                           
    877. 

  877. !                                                   units      size
    878. 

  878. !  number of soundings in zonal direction (m)       n/d        1
    879. 

  879. !  number of soundings in meridional direction (n)  n/d        1
    880. 

  880. !  maximum number of soundings in                   n/d        1
    881. 

  881. !         meridional direction (ndim>=n)
    882. 

  882. !  number of atmospheric layers (np)                n/d        1
    883. 

  883. !  level pressure (pl)                              mb     m*ndim*(np+1)
    884. 

  884. !  layer temperature (ta)                           k        m*ndim*np
    885. 

  885. !  layer specific humidity (wa)                     gm/gm    m*ndim*np
    886. 

  886. !  layer ozone concentration (oa)                   gm/gm    m*ndim*np
    887. 

  887. !  co2 mixing ratio by volumn (co2)                 pppv       1
    888. 

  888. !  option for scaling cloud optical thickness       n/d        1
    889. 

  889. !        overcast="true" if scaling is NOT required
    890. 

  890. !        overcast="fasle" if scaling is required
    891. 

  891. !  option for cloud optical thickness               n/d        1
    892. 

  892. !        cldwater="true" if cwc is provided
    893. 

  893. !        cldwater="false" if taucld is provided
    894. 

  894. !  cloud water mixing ratio (cwc)                  gm/gm     m*ndim*np*2
    895. 

  895. !        index 1 for ice particles
    896. 

  896. !        index 2 for liquid drops
    897. 

  897. !  cloud optical thickness (taucld)                 n/d      m*ndim*np*2
    898. 

  898. !        index 1 for ice particles
    899. 

  899. !        index 2 for liquid drops
    900. 

  900. !  effective cloud-particle size (reff)          micrometer m*ndim*np*2
    901. 

  901. !        index 1 for ice particles
    902. 

  902. !        index 2 for liquid drops
    903. 

  903. !  cloud amount (fcld)                            fraction   m*ndim*np
    904. 

  904. !  level index separating high and middle           n/d        1
    905. 

  905. !        clouds (ict)
    906. 

  906. !  level index separating middle and low            n/d        1
    907. 

  907. !          clouds (icb)
    908. 

  908. !  aerosol optical thickness (taual)                n/d    m*ndim*np*11
    909. 

  909. !  aerosol single-scattering albedo (ssaal)         n/d    m*ndim*np*11
    910. 

  910. !  aerosol asymmetry factor (asyal)                 n/d    m*ndim*np*11
    911. 

  911. !        in the uv region :
    912. 

  912. !           index  1 for the 0.175-0.225 micron band
    913. 

  913. !           index  2 for the 0.225-0.245; 0.260-0.280 micron band
    914. 

  914. !           index  3 for the 0.245-0.260 micron band
    915. 

  915. !           index  4 for the 0.280-0.295 micron band
    916. 

  916. !           index  5 for the 0.295-0.310 micron band
    917. 

  917. !           index  6 for the 0.310-0.320 micron band
    918. 

  918. !           index  7 for the 0.325-0.400 micron band
    919. 

  919. !        in the par region :
    920. 

  920. !           index  8 for the 0.400-0.700 micron band
    921. 

  921. !        in the infrared region :
    922. 

  922. !           index  9 for the 0.700-1.220 micron band
    923. 

  923. !           index 10 for the 1.220-2.270 micron band
    924. 

  924. !           index 11 for the 2.270-10.00 micron band
    925. 

  925. !   cosine of solar zenith angle (cosz)              n/d      m*ndim
    926. 

  926. !   uv+visible sfc albedo for beam radiation
    927. 

  927. !        for wavelengths<0.7 micron (rsuvbm)    fraction   m*ndim
    928. 

  928. !   uv+visible sfc albedo for diffuse radiation
    929. 

  929. !        for wavelengths<0.7 micron (rsuvdf)    fraction   m*ndim
    930. 

  930. !   ir sfc albedo for beam radiation
    931. 

  931. !        for wavelengths>0.7 micron  (rsirbm)   fraction   m*ndim
    932. 

  932. !   ir sfc albedo for diffuse radiation (rsirdf)   fraction   m*ndim
    933. 

  933. !
    934. 

  934. !----- Output parameters
    935. 

  935. !
    936. 

  936. !   all-sky flux (downward minus upward) (flx)    fraction m*ndim*(np+1)
    937. 

  937. !   clear-sky flux (downward minus upward) (flc)  fraction m*ndim*(np+1)
    938. 

  938. !   all-sky direct downward uv (0.175-0.4 micron)
    939. 

  939. !                flux at the surface (fdiruv)      fraction   m*ndim
    940. 

  940. !   all-sky diffuse downward uv flux at
    941. 

  941. !                the surface (fdifuv)              fraction   m*ndim
    942. 

  942. !   all-sky direct downward par (0.4-0.7 micron)
    943. 

  943. !                flux at the surface (fdirpar)     fraction   m*ndim
    944. 

  944. !   all-sky diffuse downward par flux at
    945. 

  945. !                the surface (fdifpar)             fraction   m*ndim
    946. 

  946. !   all-sky direct downward ir (0.7-10 micron)
    947. 

  947. !                flux at the surface (fdirir)      fraction   m*ndim
    948. 

  948. !   all-sky diffuse downward ir flux at
    949. 

  949. !                the surface (fdifir)              fraction   m*ndim
    950. 

  950. !
    951. 

  951. !----- Notes:
    952. 

  952. !
    953. 

  953. !    (1) The unit of "flux" is fraction of the incoming solar radiation
    954. 

  954. !        at the top of the atmosphere.  Therefore, fluxes should
    955. 

  955. !        be equal to "flux" multiplied by the extra-terrestrial solar
    956. 

  956. !        flux and the cosine of solar zenith angle.
    957. 

  957. !    (2) pl(i,j,1) is the pressure at the top of the model, and
    958. 

  958. !        pl(i,j,np+1) is the surface pressure.
    959. 

  959. !    (3) the pressure levels ict and icb correspond approximately
    960. 

  960. !        to 400 and 700 mb.
    961. 

  961. !    (4) if overcast='true', the clear-sky flux, flc, is not computed.
    962. 

  962. !        
    963. 

  963. !**************************************************************************
    964. 

  964.       implicit none
    965. 

  965. !**************************************************************************
    966. 

  966. 

  967. !-----input parameters
    968. 

  968. 

  969.       integer m,n,ndim,np
    970. 

  970.       integer ict(m,ndim),icb(m,ndim)
    971. 

  971.       real pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np)
    972. 

  972.       real cwc(m,ndim,np,2),taucld(m,ndim,np,2),reff(m,ndim,np,2), &
    973. 

  973.              fcld(m,ndim,np)
    974. 

  974.       real taual(m,ndim,np,11),ssaal(m,ndim,np,11),asyal(m,ndim,np,11)
    975. 

  975.       real cosz(m,ndim),rsuvbm(m,ndim),rsuvdf(m,ndim), &
    976. 

  976.              rsirbm(m,ndim),rsirdf(m,ndim)           
    977. 

  977.       logical overcast,cldwater
    978. 

  978. 

  979. !-----output parameters
    980. 

  980. 

  981.       real flx(m,ndim,np+1),flc(m,ndim,np+1)
    982. 

  982.       real flxu(m,ndim,np+1),flxd(m,ndim,np+1)
    983. 

  983.       real fdiruv (m,ndim),fdifuv (m,ndim)
    984. 

  984.       real fdirpar(m,ndim),fdifpar(m,ndim)
    985. 

  985.       real fdirir (m,ndim),fdifir (m,ndim)
    986. 

  986. 

  987. !-----temporary array
    988. 

  988.  
    989. 

  989.       integer i,j,k
    990. 

  990.       real cwp(m,n,np,2)
    991. 

  991.       real dp(m,n,np),wh(m,n,np),oh(m,n,np),scal(m,n,np)
    992. 

  992.       real swh(m,n,np+1),so2(m,n,np+1),df(m,n,np+1)
    993. 

  993.       real sdf(m,n),sclr(m,n),csm(m,n),x
    994. 

  994.  
    995. 

  995.       do j= 1, n 
    996. 

  996.        do i= 1, m 
    997. 

  997.           if (pl(i,j,1) .eq. 0.0) then
    998. 

  998.               pl(i,j,1)=1.0e-4
    999. 

  999.           endif
    1000. 

  1000.        enddo
    1001. 

  1001.       enddo
    1002. 

  1002. 

  1003.       do j= 1, n 
    1004. 

  1004.        do i= 1, m 
    1005. 

  1005. 

  1006.          swh(i,j,1)=0. 
    1007. 

  1007.          so2(i,j,1)=0. 
    1008. 

  1008. 

  1009. !-----csm is the effective secant of the solar zenith angle
    1010. 

  1010. !     see equation (12) of Lacis and Hansen (1974, JAS)    
    1011. 

  1011.  
    1012. 

  1012.          csm(i,j)=35./sqrt(1224.*cosz(i,j)*cosz(i,j)+1.)
    1013. 

  1013. 

  1014.        enddo 
    1015. 

  1015.       enddo
    1016. 

  1016. 

  1017.       do k= 1, np
    1018. 

  1018.        do j= 1, n
    1019. 

  1019.          do i= 1, m
    1020. 

  1020. 

  1021. !-----compute layer thickness and pressure-scaling function. 
    1022. 

  1022. !     indices for the surface level and surface layer
    1023. 

  1023. !     are np+1 and np, respectively.
    1024. 

  1024.  
    1025. 

  1025.           dp(i,j,k)=pl(i,j,k+1)-pl(i,j,k)
    1026. 

  1026.           scal(i,j,k)=dp(i,j,k)*(.5*(pl(i,j,k)+pl(i,j,k+1))/300.)**.8
    1027. 

  1027.  
    1028. 

  1028. !-----compute scaled water vapor amount, unit is g/cm**2
    1029. 

  1029. !     note: the sign prior to the constant 0.00135 was incorrectly 
    1030. 

  1030. !           set to negative in the previous version
    1031. 

  1031. 

  1032.           wh(i,j,k)=1.02*wa(i,j,k)*scal(i,j,k)* &
    1033. 

  1033.                     (1.+0.00135*(ta(i,j,k)-240.)) +1.e-11
    1034. 

  1034.           swh(i,j,k+1)=swh(i,j,k)+wh(i,j,k)
    1035. 

  1035. 

  1036. !-----compute ozone amount, unit is (cm-atm)stp
    1037. 

  1037. !     the number 466.7 is a conversion factor from g/cm**2 to (cm-atm)stp
    1038. 

  1038.  
    1039. 

  1039.           oh(i,j,k)=1.02*oa(i,j,k)*dp(i,j,k)*466.7 +1.e-11
    1040. 

  1040. 

  1041. !-----compute layer cloud water amount (gm/m**2)
    1042. 

  1042. !     the index is 1 for ice crystals and 2 for liquid drops
    1043. 

  1043. 

  1044.           cwp(i,j,k,1)=1.02*10000.*cwc(i,j,k,1)*dp(i,j,k)
    1045. 

  1045.           cwp(i,j,k,2)=1.02*10000.*cwc(i,j,k,2)*dp(i,j,k)
    1046. 

  1046. 

  1047.         enddo
    1048. 

  1048.        enddo
    1049. 

  1049.       enddo
    1050. 

  1050. 

  1051. !-----initialize fluxes for all-sky (flx), clear-sky (flc), and
    1052. 

  1052. !     flux reduction (df)
    1053. 

  1053. 

  1054.       do k=1, np+1
    1055. 

  1055.        do j=1, n
    1056. 

  1056.         do i=1, m
    1057. 

  1057.           flx(i,j,k)=0.
    1058. 

  1058.           flc(i,j,k)=0.
    1059. 

  1059.           flxu(i,j,k)=0.
    1060. 

  1060.           flxd(i,j,k)=0.
    1061. 

  1061.           df(i,j,k)=0.
    1062. 

  1062.         enddo
    1063. 

  1063.        enddo
    1064. 

  1064.       enddo
    1065. 

  1065. 

  1066. !-----compute solar uv and par fluxes
    1067. 

  1067. 

  1068.       call soluv (m,n,ndim,np,oh,dp,overcast,cldwater,  &
    1069. 

  1069.                   cwp,taucld,reff,ict,icb,fcld,cosz,    &
    1070. 

  1070.                   taual,ssaal,asyal,csm,rsuvbm,rsuvdf,  &
    1071. 

  1071.                   flx,flc,flxu,flxd,fdiruv,fdifuv,fdirpar,fdifpar)
    1072. 

  1072. 

  1073. !-----compute and update solar ir fluxes
    1074. 

  1074. 

  1075.       call solir (m,n,ndim,np,wh,overcast,cldwater,     &
    1076. 

  1076.                   cwp,taucld,reff,ict,icb,fcld,cosz,    &
    1077. 

  1077.                   taual,ssaal,asyal,csm,rsirbm,rsirdf,  &
    1078. 

  1078.                   flx,flc,flxu,flxd,fdirir,fdifir)
    1079. 

  1079. 

  1080. !-----compute scaled o2 amount, unit is (cm-atm)stp.
    1081. 

  1081. 

  1082.       do k= 1, np
    1083. 

  1083.        do j= 1, n
    1084. 

  1084.         do i= 1, m
    1085. 

  1085.           so2(i,j,k+1)=so2(i,j,k)+165.22*scal(i,j,k)
    1086. 

  1086.         enddo
    1087. 

  1087.        enddo
    1088. 

  1088.       enddo
    1089. 

  1089. 

  1090. !-----compute flux reduction due to oxygen following
    1091. 

  1091. !      chou (J. climate, 1990). The fraction 0.0287 is the
    1092. 

  1092. !      extraterrestrial solar flux in the o2 bands.
    1093. 

  1093. 

  1094.        do k= 2, np+1
    1095. 

  1095.         do j= 1, n
    1096. 

  1096.          do i= 1, m
    1097. 

  1097.            x=so2(i,j,k)*csm(i,j)
    1098. 

  1098.            df(i,j,k)=df(i,j,k)+0.0287*(1.-exp(-0.00027*sqrt(x)))
    1099. 

  1099.          enddo
    1100. 

  1100.         enddo
    1101. 

  1101.        enddo          
    1102. 

  1102. 

  1103. !-----compute scaled co2 amounts. unit is (cm-atm)stp.
    1104. 

  1104. 

  1105.       do k= 1, np
    1106. 

  1106.        do j= 1, n
    1107. 

  1107.         do i= 1, m
    1108. 

  1108.          so2(i,j,k+1)=so2(i,j,k)+co2*789.*scal(i,j,k)+1.e-11
    1109. 

  1109.         enddo
    1110. 

  1110.        enddo
    1111. 

  1111.       enddo
    1112. 

  1112. 

  1113. !-----compute and update flux reduction due to co2 following
    1114. 

  1114. !     chou (J. Climate, 1990)
    1115. 

  1115. 

  1116.       call flxco2(m,n,np,so2,swh,csm,df)
    1117. 

  1117. 

  1118. !-----adjust for the effect of o2 cnd co2 on clear-sky fluxes.
    1119. 

  1119. 

  1120.       do k= 2, np+1
    1121. 

  1121.        do j= 1, n
    1122. 

  1122.         do i= 1, m
    1123. 

  1123.           flc(i,j,k)=flc(i,j,k)-df(i,j,k)
    1124. 

  1124.         enddo
    1125. 

  1125.        enddo
    1126. 

  1126.       enddo
    1127. 

  1127. 

  1128. !-----adjust for the all-sky fluxes due to o2 and co2.  It is
    1129. 

  1129. !     assumed that o2 and co2 have no effects on solar radiation
    1130. 

  1130. !     below clouds.
    1131. 

  1131. 

  1132.       do j=1,n
    1133. 

  1133.        do i=1,m
    1134. 

  1134.         sdf(i,j)=0.0
    1135. 

  1135.         sclr(i,j)=1.0
    1136. 

  1136.        enddo
    1137. 

  1137.       enddo
    1138. 

  1138. 

  1139.       do k=1,np
    1140. 

  1140.        do j=1,n
    1141. 

  1141.         do i=1,m
    1142. 

  1142. 

  1143. !-----sclr is the fraction of clear sky.
    1144. 

  1144. !     sdf is the flux reduction below clouds.
    1145. 

  1145. 

  1146.          if(fcld(i,j,k).gt.0.01) then
    1147. 

  1147.           sdf(i,j)=sdf(i,j)+df(i,j,k)*sclr(i,j)*fcld(i,j,k)
    1148. 

  1148.           sclr(i,j)=sclr(i,j)*(1.-fcld(i,j,k))
    1149. 

  1149.          endif
    1150. 

  1150.           flx(i,j,k+1)=flx(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j)
    1151. 

  1151.           flxu(i,j,k+1)=flxu(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j)
    1152. 

  1152.           flxd(i,j,k+1)=flxd(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j) ! SG: same as flux????
    1153. 

  1153. 

  1154.         enddo
    1155. 

  1155.        enddo
    1156. 

  1156.       enddo
    1157. 

  1157. 

  1158. !-----adjustment for the direct downward ir flux.
    1159. 

  1159. 

  1160.       do j= 1, n
    1161. 

  1161.        do i= 1, m
    1162. 

  1162.         flc(i,j,np+1)=flc(i,j,np+1)+df(i,j,np+1)*rsirbm(i,j)
    1163. 

  1163.         flx(i,j,np+1)=flx(i,j,np+1)+(sdf(i,j)+ &
    1164. 

  1164.                          df(i,j,np+1)*sclr(i,j))*rsirbm(i,j)
    1165. 

  1165.         flxu(i,j,np+1)=flxu(i,j,np+1)+(sdf(i,j)+ &
    1166. 

  1166.                          df(i,j,np+1)*sclr(i,j))*rsirbm(i,j)
    1167. 

  1167.         flxd(i,j,np+1)=flxd(i,j,np+1)+(sdf(i,j)+ &
    1168. 

  1168.                          df(i,j,np+1)*sclr(i,j))*rsirbm(i,j)
    1169. 

  1169.         fdirir(i,j)=fdirir(i,j)-(sdf(i,j)+df(i,j,np+1)*sclr(i,j))
    1170. 

  1170.        enddo
    1171. 

  1171.       enddo
    1172. 

  1172. 

  1173.       end subroutine sorad 
    1174. 

  1174. 

  1175. !************************************************************************
    1176. 

  1176. 

  1177.       subroutine soluv (m,n,ndim,np,oh,dp,overcast,cldwater,            &
    1178. 

  1178.                 cwp,taucld,reff,ict,icb,fcld,cosz,                      &
    1179. 

  1179.                 taual,ssaal,asyal,csm,rsuvbm,rsuvdf,                    &
    1180. 

  1180.                 flx,flc,flxu,flxd,fdiruv,fdifuv,fdirpar,fdifpar)
    1181. 

  1181. 

  1182. !************************************************************************
    1183. 

  1183. !  compute solar fluxes in the uv+par region. the spectrum is
    1184. 

  1184. !  grouped into 8 bands:
    1185. 

  1185. !  
    1186. 

  1186. !              Band     Micrometer
    1187. 

  1187. !
    1188. 

  1188. !       UV-C    1.     .175 - .225
    1189. 

  1189. !               2.     .225 - .245
    1190. 

  1190. !                      .260 - .280
    1191. 

  1191. !               3.     .245 - .260
    1192. 

  1192. !
    1193. 

  1193. !       UV-B    4.     .280 - .295
    1194. 

  1194. !               5.     .295 - .310
    1195. 

  1195. !               6.     .310 - .320
    1196. 

  1196. !      
    1197. 

  1197. !       UV-A    7.     .320 - .400
    1198. 

  1198. !      
    1199. 

  1199. !       PAR     8.     .400 - .700
    1200. 

  1200. !
    1201. 

  1201. !----- Input parameters:                            units      size
    1202. 

  1202. !
    1203. 

  1203. !  number of soundings in zonal direction (m)       n/d        1
    1204. 

  1204. !  number of soundings in meridional direction (n)  n/d        1
    1205. 

  1205. !  maximum number of soundings in                   n/d        1
    1206. 

  1206. !        meridional direction (ndim)
    1207. 

  1207. !  number of atmospheric layers (np)                n/d        1
    1208. 

  1208. !  layer ozone content (oh)                      (cm-atm)stp m*n*np
    1209. 

  1209. !  layer pressure thickness (dp)                    mb       m*n*np
    1210. 

  1210. !  option for scaling cloud optical thickness       n/d        1
    1211. 

  1211. !        overcast="true" if scaling is NOT required
    1212. 

  1212. !        overcast="fasle" if scaling is required
    1213. 

  1213. !  input option for cloud optical thickness         n/d        1
    1214. 

  1214. !        cldwater="true" if taucld is provided
    1215. 

  1215. !        cldwater="false" if cwp is provided
    1216. 

  1216. !  cloud water amount (cwp)                        gm/m**2   m*n*np*2
    1217. 

  1217. !        index 1 for ice particles
    1218. 

  1218. !        index 2 for liquid drops
    1219. 

  1219. !  cloud optical thickness (taucld)                 n/d     m*ndim*np*2
    1220. 

  1220. !       index 1 for ice paticles
    1221. 

  1221. !       index 2 for liquid particles
    1222. 

  1222. !  effective cloud-particle size (reff)          micrometer m*ndim*np*2
    1223. 

  1223. !       index 1 for ice paticles
    1224. 

  1224. !       index 2 for liquid particles
    1225. 

  1225. !  level indiex separating high and                 n/d      m*n
    1226. 

  1226. !       middle clouds (ict)
    1227. 

  1227. !  level indiex separating middle and               n/d      m*n
    1228. 

  1228. !       low clouds (icb)
    1229. 

  1229. !  cloud amount (fcld)                            fraction   m*ndim*np
    1230. 

  1230. !  cosine of solar zenith angle (cosz)              n/d      m*ndim
    1231. 

  1231. !  aerosol optical thickness (taual)                n/d    m*ndim*np*11
    1232. 

  1232. !  aerosol single-scattering albedo (ssaal)         n/d    m*ndim*np*11
    1233. 

  1233. !  aerosol asymmetry factor (asyal)                 n/d    m*ndim*np*11
    1234. 

  1234. !  cosecant of the solar zenith angle (csm)         n/d      m*n
    1235. 

  1235. !  uv+par surface albedo for beam                 fraction   m*ndim
    1236. 

  1236. !       radiation (rsuvbm)
    1237. 

  1237. !  uv+par surface albedo for diffuse              fraction   m*ndim
    1238. 

  1238. !       radiation (rsuvdf)
    1239. 

  1239. !
    1240. 

  1240. !---- temporary array
    1241. 

  1241. !
    1242. 

  1242. !  scaled cloud optical thickness                   n/d      m*n*np
    1243. 

  1243. !       for beam radiation (tauclb)
    1244. 

  1244. !  scaled cloud optical thickness                   n/d      m*n*np
    1245. 

  1245. !       for diffuse radiation  (tauclf)     
    1246. 

  1246. !
    1247. 

  1247. !----- output (updated) parameters:
    1248. 

  1248. !
    1249. 

  1249. !  all-sky net downward flux (flx)               fraction  m*ndim*(np+1)
    1250. 

  1250. !  clear-sky net downward flux (flc)             fraction  m*ndim*(np+1)
    1251. 

  1251. !  all-sky direct downward uv flux at
    1252. 

  1252. !       the surface (fdiruv)                     fraction    m*ndim
    1253. 

  1253. !  all-sky diffuse downward uv flux at
    1254. 

  1254. !       the surface (fdifuv)                     fraction    m*ndim
    1255. 

  1255. !  all-sky direct downward par flux at
    1256. 

  1256. !       the surface (fdirpar)                    fraction    m*ndim
    1257. 

  1257. !  all-sky diffuse downward par flux at
    1258. 

  1258. !       the surface (fdifpar)                    fraction    m*ndim
    1259. 

  1259. !
    1260. 

  1260. !***********************************************************************
    1261. 

  1261.       implicit none
    1262. 

  1262. !***********************************************************************
    1263. 

  1263. 

  1264. !-----input parameters
    1265. 

  1265. 

  1266.       integer m,n,ndim,np
    1267. 

  1267.       integer ict(m,ndim),icb(m,ndim)
    1268. 

  1268.       real taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np)
    1269. 

  1269.       real cc(m,n,3),cosz(m,ndim)
    1270. 

  1270.       real cwp(m,n,np,2),oh(m,n,np),dp(m,n,np)
    1271. 

  1271.       real taual(m,ndim,np,11),ssaal(m,ndim,np,11),asyal(m,ndim,np,11)
    1272. 

  1272.       real rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n)
    1273. 

  1273.       logical overcast,cldwater
    1274. 

  1274. 

  1275. !-----output (updated) parameter
    1276. 

  1276. 

  1277.       real flx(m,ndim,np+1),flc(m,ndim,np+1)
    1278. 

  1278.       real flxu(m,ndim,np+1),flxd(m,ndim,np+1)
    1279. 

  1279.       real fdiruv (m,ndim),fdifuv (m,ndim)
    1280. 

  1280.       real fdirpar(m,ndim),fdifpar(m,ndim)
    1281. 

  1281. 

  1282. !-----static parameters
    1283. 

  1283. 

  1284.       integer nband
    1285. 

  1285.       parameter (nband=8)
    1286. 

  1286.       real hk(nband),xk(nband),ry(nband)
    1287. 

  1287.       real aig(3),awg(3)
    1288. 

  1288. 

  1289. !-----temporary array
    1290. 

  1290. 

  1291.       integer i,j,k,ib
    1292. 

  1292.       real tauclb(m,n,np),tauclf(m,n,np),asycl(m,n,np)
    1293. 

  1293.       real taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto
    1294. 

  1294.       real taux,reff1,reff2,g1,g2
    1295. 

  1295.       real td(m,n,np+1,2),rr(m,n,np+1,2),tt(m,n,np+1,2), &
    1296. 

  1296.              rs(m,n,np+1,2),ts(m,n,np+1,2)
    1297. 

  1297.       real fall(m,n,np+1),fclr(m,n,np+1),fsdir(m,n),fsdif(m,n)
    1298. 

  1298.       real fallu(m,n,np+1),falld(m,n,np+1)
    1299. 

  1299.       real asyclt(m,n)
    1300. 

  1300.       real rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
    1301. 

  1301.       real rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)
    1302. 

  1302. 

  1303. !-----hk is the fractional extra-terrestrial solar flux in each
    1304. 

  1304. !     of the 8 bands.  the sum of hk is 0.47074.
    1305. 

  1305. 

  1306.       data hk/.00057, .00367, .00083, .00417,  &
    1307. 

  1307.               .00600, .00556, .05913, .39081/
    1308. 

  1308. 

  1309. !-----xk is the ozone absorption coefficient. unit: /(cm-atm)stp
    1310. 

  1310. 

  1311.       data xk /30.47, 187.2,  301.9,   42.83, &
    1312. 

  1312.                7.09,  1.25,   0.0345,  0.0539/
    1313. 

  1313. 

  1314. !-----ry is the extinction coefficient for Rayleigh scattering.
    1315. 

  1315. !     unit: /mb.
    1316. 

  1316. 

  1317.       data ry /.00604, .00170, .00222, .00132, &
    1318. 

  1318.                .00107, .00091, .00055, .00012/
    1319. 

  1319. 

  1320. !-----coefficients for computing the asymmetry factor of ice clouds
    1321. 

  1321. !     from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2, independent
    1322. 

  1322. !     of spectral band.
    1323. 

  1323. 

  1324.       data aig/.74625000,.00105410,-.00000264/
    1325. 

  1325. 

  1326. !-----coefficients for computing the asymmetry factor of liquid
    1327. 

  1327. !     clouds from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2,
    1328. 

  1328. !     independent of spectral band.
    1329. 

  1329. 

  1330.       data awg/.82562000,.00529000,-.00014866/
    1331. 

  1331. 

  1332. !-----initialize fdiruv, fdifuv, surface reflectances and transmittances.
    1333. 

  1333. !     cc is the maximum cloud cover in each of the three cloud groups.
    1334. 

  1334.             
    1335. 

  1335.       do j= 1, n
    1336. 

  1336.        do i= 1, m                    
    1337. 

  1337.          fdiruv(i,j)=0.0
    1338. 

  1338.          fdifuv(i,j)=0.0
    1339. 

  1339.          rr(i,j,np+1,1)=rsuvbm(i,j)
    1340. 

  1340.          rr(i,j,np+1,2)=rsuvbm(i,j)
    1341. 

  1341.          rs(i,j,np+1,1)=rsuvdf(i,j)
    1342. 

  1342.          rs(i,j,np+1,2)=rsuvdf(i,j)
    1343. 

  1343.          td(i,j,np+1,1)=0.0
    1344. 

  1344.          td(i,j,np+1,2)=0.0
    1345. 

  1345.          tt(i,j,np+1,1)=0.0
    1346. 

  1346.          tt(i,j,np+1,2)=0.0
    1347. 

  1347.          ts(i,j,np+1,1)=0.0
    1348. 

  1348.          ts(i,j,np+1,2)=0.0
    1349. 

  1349.          cc(i,j,1)=0.0
    1350. 

  1350.          cc(i,j,2)=0.0
    1351. 

  1351.          cc(i,j,3)=0.0
    1352. 

  1352.        enddo
    1353. 

  1353.       enddo
    1354. 

  1354. 

  1355. 

  1356. !-----compute cloud optical thickness
    1357. 

  1357. 

  1358.       if (cldwater) then
    1359. 

  1359. 

  1360.        do k= 1, np
    1361. 

  1361.         do j= 1, n
    1362. 

  1362.          do i= 1, m
    1363. 

  1363.           taucld(i,j,k,1)=cwp(i,j,k,1)*( 3.33e-4+2.52/reff(i,j,k,1))
    1364. 

  1364.           taucld(i,j,k,2)=cwp(i,j,k,2)*(-6.59e-3+1.65/reff(i,j,k,2))
    1365. 

  1365.          enddo
    1366. 

  1366.         enddo
    1367. 

  1367.        enddo
    1368. 

  1368. 

  1369.       endif
    1370. 

  1370. 

  1371. !-----options for scaling cloud optical thickness
    1372. 

  1372. 

  1373.       if (overcast) then
    1374. 

  1374. 

  1375.        do k= 1, np
    1376. 

  1376.         do j= 1, n
    1377. 

  1377.          do i= 1, m
    1378. 

  1378.           tauclb(i,j,k)=taucld(i,j,k,1)+taucld(i,j,k,2)
    1379. 

  1379.           tauclf(i,j,k)=tauclb(i,j,k)
    1380. 

  1380.          enddo
    1381. 

  1381.         enddo
    1382. 

  1382.        enddo
    1383. 

  1383. 

  1384.        do k= 1, 3
    1385. 

  1385.         do j= 1, n
    1386. 

  1386.          do i= 1, m
    1387. 

  1387.            cc(i,j,k)=1.0
    1388. 

  1388.          enddo
    1389. 

  1389.         enddo
    1390. 

  1390.        enddo
    1391. 

  1391. 

  1392.       else
    1393. 

  1393. 

  1394. !-----scale cloud optical thickness in each layer from taucld (with
    1395. 

  1395. !     cloud amount fcld) to tauclb and tauclf (with cloud amount cc).
    1396. 

  1396. !     tauclb is the scaled optical thickness for beam radiation and
    1397. 

  1397. !     tauclf is for diffuse radiation.
    1398. 

  1398. 

  1399.        call cldscale(m,n,ndim,np,cosz,fcld,taucld,ict,icb,  &
    1400. 

  1400.                     cc,tauclb,tauclf)
    1401. 

  1401. 

  1402.       endif
    1403. 

  1403. 

  1404. !-----compute cloud asymmetry factor for a mixture of
    1405. 

  1405. !     liquid and ice particles.  unit of reff is micrometers.
    1406. 

  1406. 

  1407.       do k= 1, np
    1408. 

  1408. 

  1409.        do j= 1, n
    1410. 

  1410.         do i= 1, m
    1411. 

  1411. 

  1412.            asyclt(i,j)=1.0
    1413. 

  1413. 

  1414.            taux=taucld(i,j,k,1)+taucld(i,j,k,2)
    1415. 

  1415.           if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
    1416. 

  1416. 

  1417.            reff1=min(reff(i,j,k,1),130.)
    1418. 

  1418.            reff2=min(reff(i,j,k,2),20.0)
    1419. 

  1419. 

  1420.            g1=(aig(1)+(aig(2)+aig(3)*reff1)*reff1)*taucld(i,j,k,1)
    1421. 

  1421.            g2=(awg(1)+(awg(2)+awg(3)*reff2)*reff2)*taucld(i,j,k,2)
    1422. 

  1422.            asyclt(i,j)=(g1+g2)/taux
    1423. 

  1423. 

  1424.           endif
    1425. 

  1425. 

  1426.         enddo
    1427. 

  1427.        enddo
    1428. 

  1428. 

  1429.        do j=1,n
    1430. 

  1430.         do i=1,m
    1431. 

  1431.            asycl(i,j,k)=asyclt(i,j)
    1432. 

  1432.         enddo
    1433. 

  1433.        enddo
    1434. 

  1434. 

  1435.       enddo
    1436. 

  1436. 

  1437. !-----integration over spectral bands
    1438. 

  1438. 

  1439.       do 100 ib=1,nband
    1440. 

  1440. 

  1441.        do 300 k= 1, np
    1442. 

  1442. 

  1443.         do j= 1, n
    1444. 

  1444.          do i= 1, m
    1445. 

  1445. 

  1446. !-----compute ozone and rayleigh optical thicknesses
    1447. 

  1447. 

  1448.           taurs=ry(ib)*dp(i,j,k)
    1449. 

  1449.           tauoz=xk(ib)*oh(i,j,k)
    1450. 

  1450.  
    1451. 

  1451. !-----compute clear-sky optical thickness, single scattering albedo,
    1452. 

  1452. !     and asymmetry factor
    1453. 

  1453. 

  1454.           tausto=taurs+tauoz+taual(i,j,k,ib)+1.0e-8
    1455. 

  1455.           ssatau=ssaal(i,j,k,ib)*taual(i,j,k,ib)+taurs
    1456. 

  1456.           asysto=asyal(i,j,k,ib)*ssaal(i,j,k,ib)*taual(i,j,k,ib)
    1457. 

  1457. 

  1458.           tauto=tausto
    1459. 

  1459.           ssato=ssatau/tauto+1.0e-8
    1460. 

  1460.           ssato=min(ssato,0.999999)
    1461. 

  1461.           asyto=asysto/(ssato*tauto)
    1462. 

  1462. 

  1463. !-----compute reflectance and transmittance for cloudless layers
    1464. 

  1464. 

  1465. !-                 for direct incident radiation
    1466. 

  1466. 

  1467.           call deledd (tauto,ssato,asyto,csm(i,j),  &
    1468. 

  1468.                        rr1t(i,j),tt1t(i,j),td1t(i,j))
    1469. 

  1469. 

  1470. !-                 for diffuse incident radiation
    1471. 

  1471. 

  1472.           call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j))
    1473. 

  1473. 

  1474. !-----compute reflectance and transmittance for cloud layers
    1475. 

  1475. 

  1476.          if (tauclb(i,j,k).lt.0.01 .or. fcld(i,j,k).lt.0.01) then
    1477. 

  1477. 

  1478.           rr2t(i,j)=rr1t(i,j)
    1479. 

  1479.           tt2t(i,j)=tt1t(i,j)
    1480. 

  1480.           td2t(i,j)=td1t(i,j)
    1481. 

  1481.           rs2t(i,j)=rs1t(i,j)
    1482. 

  1482.           ts2t(i,j)=ts1t(i,j)
    1483. 

  1483. 

  1484.          else
    1485. 

  1485. 

  1486. !--                for direct incident radiation
    1487. 

  1487. 

  1488.           tauto=tausto+tauclb(i,j,k)
    1489. 

  1489.           ssato=(ssatau+tauclb(i,j,k))/tauto+1.0e-8
    1490. 

  1490.           ssato=min(ssato,0.999999)
    1491. 

  1491.           asyto=(asysto+asycl(i,j,k)*tauclb(i,j,k))/(ssato*tauto)
    1492. 

  1492. 

  1493.           call deledd (tauto,ssato,asyto,csm(i,j),  &
    1494. 

  1494.                        rr2t(i,j),tt2t(i,j),td2t(i,j))
    1495. 

  1495. 

  1496. !--                for diffuse incident radiation
    1497. 

  1497. 

  1498.           tauto=tausto+tauclf(i,j,k)
    1499. 

  1499.           ssato=(ssatau+tauclf(i,j,k))/tauto+1.0e-8
    1500. 

  1500.           ssato=min(ssato,0.999999)
    1501. 

  1501.           asyto=(asysto+asycl(i,j,k)*tauclf(i,j,k))/(ssato*tauto)
    1502. 

  1502. 

  1503.           call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))
    1504. 

  1504. 

  1505.          endif
    1506. 

  1506. 

  1507.         enddo
    1508. 

  1508.        enddo
    1509. 

  1509. 

  1510.         do j=1,n
    1511. 

  1511.          do i=1,m
    1512. 

  1512.             rr(i,j,k,1)=rr1t(i,j)
    1513. 

  1513.          enddo
    1514. 

  1514.         enddo
    1515. 

  1515.         do j=1,n
    1516. 

  1516.          do i=1,m
    1517. 

  1517.             tt(i,j,k,1)=tt1t(i,j)
    1518. 

  1518.          enddo
    1519. 

  1519.         enddo
    1520. 

  1520.         do j=1,n
    1521. 

  1521.          do i=1,m
    1522. 

  1522.             td(i,j,k,1)=td1t(i,j)
    1523. 

  1523.          enddo
    1524. 

  1524.         enddo
    1525. 

  1525.         do j=1,n
    1526. 

  1526.          do i=1,m
    1527. 

  1527.             rs(i,j,k,1)=rs1t(i,j)
    1528. 

  1528.          enddo
    1529. 

  1529.         enddo
    1530. 

  1530.         do j=1,n
    1531. 

  1531.          do i=1,m
    1532. 

  1532.             ts(i,j,k,1)=ts1t(i,j)
    1533. 

  1533.          enddo
    1534. 

  1534.         enddo
    1535. 

  1535. 

  1536.         do j=1,n
    1537. 

  1537.          do i=1,m
    1538. 

  1538.             rr(i,j,k,2)=rr2t(i,j)
    1539. 

  1539.          enddo
    1540. 

  1540.         enddo
    1541. 

  1541.         do j=1,n
    1542. 

  1542.          do i=1,m
    1543. 

  1543.             tt(i,j,k,2)=tt2t(i,j)
    1544. 

  1544.          enddo
    1545. 

  1545.         enddo
    1546. 

  1546.         do j=1,n
    1547. 

  1547.          do i=1,m
    1548. 

  1548.             td(i,j,k,2)=td2t(i,j)
    1549. 

  1549.          enddo
    1550. 

  1550.         enddo
    1551. 

  1551.         do j=1,n
    1552. 

  1552.          do i=1,m
    1553. 

  1553.             rs(i,j,k,2)=rs2t(i,j)
    1554. 

  1554.          enddo
    1555. 

  1555.         enddo
    1556. 

  1556.         do j=1,n
    1557. 

  1557.          do i=1,m
    1558. 

  1558.             ts(i,j,k,2)=ts2t(i,j)
    1559. 

  1559.          enddo
    1560. 

  1560.         enddo
    1561. 

  1561. 

  1562.  300  continue
    1563. 

  1563. 

  1564. !-----flux calculations
    1565. 

  1565.  
    1566. 

  1566.         call cldflx (m,n,np,ict,icb,overcast,cc,rr,tt,td,rs,ts, &
    1567. 

  1567.                      fclr,fall,fallu,falld,fsdir,fsdif)
    1568. 

  1568. 

  1569.        do k= 1, np+1
    1570. 

  1570.         do j= 1, n
    1571. 

  1571.          do i= 1, m
    1572. 

  1572.           flx(i,j,k)=flx(i,j,k)+fall(i,j,k)*hk(ib)
    1573. 

  1573.           flxu(i,j,k)=flxu(i,j,k)+fallu(i,j,k)*hk(ib)
    1574. 

  1574.           flxd(i,j,k)=flxd(i,j,k)+falld(i,j,k)*hk(ib)
    1575. 

  1575.          enddo
    1576. 

  1576.         enddo
    1577. 

  1577.         do j= 1, n
    1578. 

  1578.          do i= 1, m
    1579. 

  1579.           flc(i,j,k)=flc(i,j,k)+fclr(i,j,k)*hk(ib)
    1580. 

  1580.          enddo
    1581. 

  1581.         enddo
    1582. 

  1582.        enddo
    1583. 

  1583. 

  1584. !-----compute downward surface fluxes in the UV and par regions
    1585. 

  1585. 

  1586.        if(ib.lt.8) then
    1587. 

  1587.          do j=1,n
    1588. 

  1588.           do i=1,m
    1589. 

  1589.            fdiruv(i,j)=fdiruv(i,j)+fsdir(i,j)*hk(ib)
    1590. 

  1590.            fdifuv(i,j)=fdifuv(i,j)+fsdif(i,j)*hk(ib)
    1591. 

  1591.          enddo
    1592. 

  1592.         enddo
    1593. 

  1593.        else
    1594. 

  1594.          do j=1,n
    1595. 

  1595.           do i=1,m
    1596. 

  1596.            fdirpar(i,j)=fsdir(i,j)*hk(ib)
    1597. 

  1597.            fdifpar(i,j)=fsdif(i,j)*hk(ib)
    1598. 

  1598.          enddo
    1599. 

  1599.         enddo
    1600. 

  1600.        endif
    1601. 

  1601. 

  1602.  100  continue
    1603. 

  1603. 

  1604.       end subroutine soluv
    1605. 

  1605. 

  1606. !************************************************************************
    1607. 

  1607. 

  1608.       subroutine solir (m,n,ndim,np,wh,overcast,cldwater,               &
    1609. 

  1609.                         cwp,taucld,reff,ict,icb,fcld,cosz,              &
    1610. 

  1610.                         taual,ssaal,asyal,csm,rsirbm,rsirdf,            &
    1611. 

  1611.                         flx,flc,flxu,flxd,fdirir,fdifir)
    1612. 

  1612. 

  1613. !************************************************************************
    1614. 

  1614. !  compute solar flux in the infrared region. The spectrum is divided
    1615. 

  1615. !   into three bands:
    1616. 

  1616. !
    1617. 

  1617. !          band   wavenumber(/cm)  wavelength (micron)
    1618. 

  1618. !          1( 9)    14300-8200         0.70-1.22
    1619. 

  1619. !          2(10)     8200-4400         1.22-2.27
    1620. 

  1620. !          3(11)     4400-1000         2.27-10.0
    1621. 

  1621. !
    1622. 

  1622. !----- Input parameters:                            units      size
    1623. 

  1623. !
    1624. 

  1624. !  number of soundings in zonal direction (m)       n/d        1
    1625. 

  1625. !  number of soundings in meridional direction (n)  n/d        1
    1626. 

  1626. !  maximum number of soundings in                   n/d        1
    1627. 

  1627. !         meridional direction (ndim)
    1628. 

  1628. !  number of atmospheric layers (np)                n/d        1
    1629. 

  1629. !  layer scaled-water vapor content (wh)          gm/cm^2    m*n*np
    1630. 

  1630. !  option for scaling cloud optical thickness       n/d        1
    1631. 

  1631. !        overcast="true" if scaling is NOT required
    1632. 

  1632. !        overcast="fasle" if scaling is required
    1633. 

  1633. !  input option for cloud optical thickness         n/d        1
    1634. 

  1634. !        cldwater="true" if taucld is provided
    1635. 

  1635. !        cldwater="false" if cwp is provided
    1636. 

  1636. !  cloud water concentration (cwp)                gm/m**2    m*n*np*2
    1637. 

  1637. !        index 1 for ice particles
    1638. 

  1638. !        index 2 for liquid drops
    1639. 

  1639. !  cloud optical thickness (taucld)                 n/d      m*ndim*np*2
    1640. 

  1640. !        index 1 for ice paticles
    1641. 

  1641. !  effective cloud-particle size (reff)           micrometer m*ndim*np*2
    1642. 

  1642. !        index 1 for ice paticles
    1643. 

  1643. !        index 2 for liquid particles
    1644. 

  1644. !  level index separating high and                  n/d      m*n
    1645. 

  1645. !        middle clouds (ict)
    1646. 

  1646. !  level index separating middle and                n/d      m*n
    1647. 

  1647. !        low clouds (icb)
    1648. 

  1648. !  cloud amount (fcld)                            fraction   m*ndim*np
    1649. 

  1649. !  aerosol optical thickness (taual)                n/d      m*ndim*np*11
    1650. 

  1650. !  aerosol single-scattering albedo (ssaal)         n/d      m*ndim*np*11
    1651. 

  1651. !  aerosol asymmetry factor (asyal)                 n/d      m*ndim*np*11 
    1652. 

  1652. !  cosecant of the solar zenith angle (csm)         n/d      m*n
    1653. 

  1653. !  near ir surface albedo for beam                fraction   m*ndim
    1654. 

  1654. !        radiation (rsirbm)
    1655. 

  1655. !  near ir surface albedo for diffuse             fraction   m*ndim
    1656. 

  1656. !        radiation (rsirdf)
    1657. 

  1657. !
    1658. 

  1658. !---- temporary array
    1659. 

  1659. !
    1660. 

  1660. !  scaled cloud optical thickness                   n/d      m*n*np
    1661. 

  1661. !          for beam radiation (tauclb)
    1662. 

  1662. !  scaled cloud optical thickness                   n/d      m*n*np
    1663. 

  1663. !          for diffuse radiation  (tauclf)     
    1664. 

  1664. !
    1665. 

  1665. !----- output (updated) parameters:
    1666. 

  1666. !
    1667. 

  1667. !  all-sky flux (downward-upward) (flx)           fraction   m*ndim*(np+1)
    1668. 

  1668. !  clear-sky flux (downward-upward) (flc)         fraction   m*ndim*(np+1)
    1669. 

  1669. !  all-sky direct downward ir flux at
    1670. 

  1670. !          the surface (fdirir)                    fraction   m*ndim
    1671. 

  1671. !  all-sky diffuse downward ir flux at
    1672. 

  1672. !          the surface (fdifir)                    fraction   m*ndim
    1673. 

  1673. !
    1674. 

  1674. !**********************************************************************
    1675. 

  1675.       implicit none
    1676. 

  1676. !**********************************************************************
    1677. 

  1677. 

  1678. !-----input parameters
    1679. 

  1679. 

  1680.       integer m,n,ndim,np
    1681. 

  1681.       integer ict(m,ndim),icb(m,ndim)
    1682. 

  1682.       real cwp(m,n,np,2),taucld(m,ndim,np,2),reff(m,ndim,np,2)
    1683. 

  1683.       real fcld(m,ndim,np),cc(m,n,3),cosz(m,ndim)
    1684. 

  1684.       real rsirbm(m,ndim),rsirdf(m,ndim)
    1685. 

  1685.       real taual(m,ndim,np,11),ssaal(m,ndim,np,11),asyal(m,ndim,np,11)
    1686. 

  1686.       real wh(m,n,np),csm(m,n)
    1687. 

  1687.       logical overcast,cldwater
    1688. 

  1688. 

  1689. !-----output (updated) parameters
    1690. 

  1690. 

  1691.       real flx(m,ndim,np+1),flc(m,ndim,np+1)
    1692. 

  1692.       real flxu(m,ndim,np+1),flxd(m,ndim,np+1)
    1693. 

  1693.       real fdirir(m,ndim),fdifir(m,ndim)
    1694. 

  1694. 

  1695. !-----static parameters
    1696. 

  1696. 

  1697.       integer nk,nband
    1698. 

  1698.       parameter (nk=10,nband=3)
    1699. 

  1699.       real xk(nk),hk(nband,nk),aib(nband,2),awb(nband,2)
    1700. 

  1700.       real aia(nband,3),awa(nband,3),aig(nband,3),awg(nband,3)
    1701. 

  1701. 

  1702. !-----temporary array
    1703. 

  1703. 

  1704.       integer ib,iv,ik,i,j,k
    1705. 

  1705.       real tauclb(m,n,np),tauclf(m,n,np)
    1706. 

  1706.       real ssacl(m,n,np),asycl(m,n,np)
    1707. 

  1707.       real rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2), &
    1708. 

  1708.              rs(m,n,np+1,2),ts(m,n,np+1,2)
    1709. 

  1709.       real fall(m,n,np+1),fclr(m,n,np+1)
    1710. 

  1710.       real fallu(m,n,np+1),falld(m,n,np+1)
    1711. 

  1711.       real fsdir(m,n),fsdif(m,n)
    1712. 

  1712. 

  1713.       real tauwv,tausto,ssatau,asysto,tauto,ssato,asyto
    1714. 

  1714.       real taux,reff1,reff2,w1,w2,g1,g2
    1715. 

  1715.       real ssaclt(m,n),asyclt(m,n)
    1716. 

  1716.       real rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n)
    1717. 

  1717.       real rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n)
    1718. 

  1718. 

  1719. !-----water vapor absorption coefficient for 10 k-intervals.
    1720. 

  1720. !     unit: cm^2/gm
    1721. 

  1721. 

  1722.       data xk/  &
    1723. 

  1723.         0.0010, 0.0133, 0.0422, 0.1334, 0.4217, &
    1724. 

  1724.         1.334,  5.623,  31.62,  177.8,  1000.0/  
    1725. 

  1725. 

  1726. !-----water vapor k-distribution function,
    1727. 

  1727. !     the sum of hk is 0.52926. unit: fraction
    1728. 

  1728. 

  1729.       data hk/  &
    1730. 

  1730.        .20673,.08236,.01074,  .03497,.01157,.00360, &
    1731. 

  1731.        .03011,.01133,.00411,  .02260,.01143,.00421, &
    1732. 

  1732.        .01336,.01240,.00389,  .00696,.01258,.00326, &
    1733. 

  1733.        .00441,.01381,.00499,  .00115,.00650,.00465, &
    1734. 

  1734.        .00026,.00244,.00245,  .00000,.00094,.00145/
    1735. 

  1735. 

  1736. !-----coefficients for computing the extinction coefficient of
    1737. 

  1737. !     ice clouds from b=aib(*,1)+aib(*,2)/reff
    1738. 

  1738. 

  1739.       data aib/ &
    1740. 

  1740.         .000333, .000333, .000333, &
    1741. 

  1741.            2.52,    2.52,    2.52/
    1742. 

  1742. 

  1743. !-----coefficients for computing the extinction coefficient of
    1744. 

  1744. !     water clouds from b=awb(*,1)+awb(*,2)/reff
    1745. 

  1745. 

  1746.       data awb/ &
    1747. 

  1747.         -0.0101, -0.0166, -0.0339, &
    1748. 

  1748.            1.72,    1.85,    2.16/
    1749. 

  1749. 

  1750. 

  1751. !-----coefficients for computing the single scattering albedo of
    1752. 

  1752. !     ice clouds from ssa=1-(aia(*,1)+aia(*,2)*reff+aia(*,3)*reff**2)
    1753. 

  1753. 

  1754.       data aia/ &
    1755. 

  1755.        -.00000260, .00215346, .08938331, &
    1756. 

  1756.         .00000746, .00073709, .00299387, &
    1757. 

  1757.         .00000000,-.00000134,-.00001038/
    1758. 

  1758. 

  1759. !-----coefficients for computing the single scattering albedo of
    1760. 

  1760. !     liquid clouds from ssa=1-(awa(*,1)+awa(*,2)*reff+awa(*,3)*reff**2)
    1761. 

  1761. 

  1762.       data awa/ &
    1763. 

  1763.         .00000007,-.00019934, .01209318, &
    1764. 

  1764.         .00000845, .00088757, .01784739, &
    1765. 

  1765.        -.00000004,-.00000650,-.00036910/
    1766. 

  1766. 

  1767. !-----coefficients for computing the asymmetry factor of ice clouds
    1768. 

  1768. !     from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2
    1769. 

  1769. 

  1770.       data aig/ &
    1771. 

  1771.         .74935228, .76098937, .84090400, &
    1772. 

  1772.         .00119715, .00141864, .00126222, &
    1773. 

  1773.        -.00000367,-.00000396,-.00000385/
    1774. 

  1774. 

  1775. !-----coefficients for computing the asymmetry factor of liquid clouds
    1776. 

  1776. !     from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2
    1777. 

  1777. 

  1778.       data awg/ &
    1779. 

  1779.         .79375035, .74513197, .83530748, &
    1780. 

  1780.         .00832441, .01370071, .00257181, &
    1781. 

  1781.        -.00023263,-.00038203, .00005519/
    1782. 

  1782. 

  1783. !-----initialize surface fluxes, reflectances, and transmittances.
    1784. 

  1784. !     cc is the maximum cloud cover in each of the three cloud groups.
    1785. 

  1785. 

  1786.       do j= 1, n
    1787. 

  1787.        do i= 1, m
    1788. 

  1788.          fdirir(i,j)=0.0
    1789. 

  1789.          fdifir(i,j)=0.0
    1790. 

  1790.          rr(i,j,np+1,1)=rsirbm(i,j)
    1791. 

  1791.          rr(i,j,np+1,2)=rsirbm(i,j)
    1792. 

  1792.          rs(i,j,np+1,1)=rsirdf(i,j)
    1793. 

  1793.          rs(i,j,np+1,2)=rsirdf(i,j)
    1794. 

  1794.          td(i,j,np+1,1)=0.0
    1795. 

  1795.          td(i,j,np+1,2)=0.0
    1796. 

  1796.          tt(i,j,np+1,1)=0.0
    1797. 

  1797.          tt(i,j,np+1,2)=0.0
    1798. 

  1798.          ts(i,j,np+1,1)=0.0
    1799. 

  1799.          ts(i,j,np+1,2)=0.0
    1800. 

  1800.          cc(i,j,1)=0.0
    1801. 

  1801.          cc(i,j,2)=0.0
    1802. 

  1802.          cc(i,j,3)=0.0
    1803. 

  1803.        enddo
    1804. 

  1804.       enddo
    1805. 

  1805. 

  1806. !-----integration over spectral bands
    1807. 

  1807. 

  1808.       do 100 ib=1,nband
    1809. 

  1809. 

  1810.        iv=ib+8
    1811. 

  1811. 

  1812. !-----compute cloud optical thickness
    1813. 

  1813. 

  1814.       if (cldwater) then
    1815. 

  1815. 

  1816.        do k= 1, np
    1817. 

  1817.         do j= 1, n
    1818. 

  1818.          do i= 1, m
    1819. 

  1819.           taucld(i,j,k,1)=cwp(i,j,k,1)*(aib(ib,1) &
    1820. 

  1820.                           +aib(ib,2)/reff(i,j,k,1))
    1821. 

  1821.           taucld(i,j,k,2)=cwp(i,j,k,2)*(awb(ib,1) &
    1822. 

  1822.                           +awb(ib,2)/reff(i,j,k,2))
    1823. 

  1823.          enddo
    1824. 

  1824.         enddo
    1825. 

  1825.        enddo
    1826. 

  1826. 

  1827.       endif
    1828. 

  1828. 

  1829. !-----options for scaling cloud optical thickness
    1830. 

  1830. 

  1831.       if (overcast) then
    1832. 

  1832. 

  1833.        do k= 1, np
    1834. 

  1834.         do j= 1, n
    1835. 

  1835.          do i= 1, m
    1836. 

  1836.           tauclb(i,j,k)=taucld(i,j,k,1)+taucld(i,j,k,2)
    1837. 

  1837.           tauclf(i,j,k)=tauclb(i,j,k)
    1838. 

  1838.          enddo
    1839. 

  1839.         enddo
    1840. 

  1840.        enddo
    1841. 

  1841. 

  1842.        do k= 1, 3
    1843. 

  1843.         do j= 1, n
    1844. 

  1844.          do i= 1, m
    1845. 

  1845.            cc(i,j,k)=1.0
    1846. 

  1846.          enddo
    1847. 

  1847.         enddo
    1848. 

  1848.        enddo
    1849. 

  1849. 

  1850.       else
    1851. 

  1851. 

  1852. !-----scale cloud optical thickness in each layer from taucld (with
    1853. 

  1853. !     cloud amount fcld) to tauclb and tauclf (with cloud amount cc).
    1854. 

  1854. !     tauclb is the scaled optical thickness for beam radiation and
    1855. 

  1855. !     tauclf is for diffuse radiation.
    1856. 

  1856. 

  1857.        call cldscale(m,n,ndim,np,cosz,fcld,taucld,ict,icb, &
    1858. 

  1858.                     cc,tauclb,tauclf)
    1859. 

  1859. 

  1860.       endif
    1861. 

  1861. 

  1862. !-----compute cloud single scattering albedo and asymmetry factor
    1863. 

  1863. !     for a mixture of ice and liquid particles.
    1864. 

  1864. 

  1865.        do k= 1, np
    1866. 

  1866. 

  1867.         do j= 1, n
    1868. 

  1868.          do i= 1, m
    1869. 

  1869. 

  1870.            ssaclt(i,j)=1.0
    1871. 

  1871.            asyclt(i,j)=1.0
    1872. 

  1872. 

  1873.            taux=taucld(i,j,k,1)+taucld(i,j,k,2)
    1874. 

  1874.           if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
    1875. 

  1875. 

  1876.            reff1=min(reff(i,j,k,1),130.)
    1877. 

  1877.            reff2=min(reff(i,j,k,2),20.0)
    1878. 

  1878. 

  1879.            w1=(1.-(aia(ib,1)+(aia(ib,2)+ &
    1880. 

  1880.                aia(ib,3)*reff1)*reff1))*taucld(i,j,k,1)
    1881. 

  1881.            w2=(1.-(awa(ib,1)+(awa(ib,2)+ &
    1882. 

  1882.                awa(ib,3)*reff2)*reff2))*taucld(i,j,k,2)
    1883. 

  1883.            ssaclt(i,j)=(w1+w2)/taux
    1884. 

  1884. 

  1885.            g1=(aig(ib,1)+(aig(ib,2)+aig(ib,3)*reff1)*reff1)*w1
    1886. 

  1886.            g2=(awg(ib,1)+(awg(ib,2)+awg(ib,3)*reff2)*reff2)*w2
    1887. 

  1887.            asyclt(i,j)=(g1+g2)/(w1+w2)
    1888. 

  1888. 

  1889.           endif
    1890. 

  1890. 

  1891.          enddo
    1892. 

  1892.         enddo
    1893. 

  1893. 

  1894.         do j=1,n
    1895. 

  1895.          do i=1,m
    1896. 

  1896.             ssacl(i,j,k)=ssaclt(i,j)
    1897. 

  1897.          enddo
    1898. 

  1898.         enddo
    1899. 

  1899.         do j=1,n
    1900. 

  1900.          do i=1,m
    1901. 

  1901.             asycl(i,j,k)=asyclt(i,j)
    1902. 

  1902.          enddo
    1903. 

  1903.         enddo
    1904. 

  1904. 

  1905.        enddo
    1906. 

  1906. 

  1907. !-----integration over the k-distribution function
    1908. 

  1908. 

  1909.          do 200 ik=1,nk
    1910. 

  1910. 

  1911.           do 300 k= 1, np
    1912. 

  1912. 

  1913.            do j= 1, n
    1914. 

  1914.             do i= 1, m
    1915. 

  1915. 

  1916.              tauwv=xk(ik)*wh(i,j,k)
    1917. 

  1917.  
    1918. 

  1918. !-----compute clear-sky optical thickness, single scattering albedo,
    1919. 

  1919. !     and asymmetry factor.
    1920. 

  1920.  
    1921. 

  1921.              tausto=tauwv+taual(i,j,k,iv)+1.0e-8
    1922. 

  1922.              ssatau=ssaal(i,j,k,iv)*taual(i,j,k,iv)
    1923. 

  1923.              asysto=asyal(i,j,k,iv)*ssaal(i,j,k,iv)*taual(i,j,k,iv)
    1924. 

  1924.  
    1925. 

  1925. !-----compute reflectance and transmittance for cloudless layers
    1926. 

  1926. 

  1927.              tauto=tausto
    1928. 

  1928.              ssato=ssatau/tauto+1.0e-8
    1929. 

  1929. 

  1930.             if (ssato .gt. 0.001) then
    1931. 

  1931. 

  1932.              ssato=min(ssato,0.999999)
    1933. 

  1933.              asyto=asysto/(ssato*tauto)
    1934. 

  1934. 

  1935. !-                 for direct incident radiation
    1936. 

  1936. 

  1937.              call deledd (tauto,ssato,asyto,csm(i,j),  &
    1938. 

  1938.                           rr1t(i,j),tt1t(i,j),td1t(i,j))
    1939. 

  1939. 

  1940. !-                 for diffuse incident radiation
    1941. 

  1941. 

  1942.              call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j))
    1943. 

  1943. 

  1944.             else
    1945. 

  1945. 

  1946.              td1t(i,j)=exp(-tauto*csm(i,j))
    1947. 

  1947.              ts1t(i,j)=exp(-1.66*tauto)
    1948. 

  1948.              tt1t(i,j)=0.0
    1949. 

  1949.              rr1t(i,j)=0.0
    1950. 

  1950.              rs1t(i,j)=0.0
    1951. 

  1951. 

  1952.             endif
    1953. 

  1953. 

  1954. !-----compute reflectance and transmittance for cloud layers
    1955. 

  1955. 

  1956.             if (tauclb(i,j,k).lt.0.01 .or. fcld(i,j,k).lt.0.01) then
    1957. 

  1957. 

  1958.              rr2t(i,j)=rr1t(i,j)
    1959. 

  1959.              tt2t(i,j)=tt1t(i,j)
    1960. 

  1960.              td2t(i,j)=td1t(i,j)
    1961. 

  1961.              rs2t(i,j)=rs1t(i,j)
    1962. 

  1962.              ts2t(i,j)=ts1t(i,j)
    1963. 

  1963. 

  1964.             else
    1965. 

  1965. 

  1966. !-                 for direct incident radiation
    1967. 

  1967. 

  1968.              tauto=tausto+tauclb(i,j,k)
    1969. 

  1969.              ssato=(ssatau+ssacl(i,j,k)*tauclb(i,j,k))/tauto+1.0e-8
    1970. 

  1970.              ssato=min(ssato,0.999999)
    1971. 

  1971.              asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclb(i,j,k))/ &
    1972. 

  1972.                    (ssato*tauto)
    1973. 

  1973. 

  1974.              call deledd (tauto,ssato,asyto,csm(i,j),  &
    1975. 

  1975.                           rr2t(i,j),tt2t(i,j),td2t(i,j))
    1976. 

  1976. 

  1977. !-                 for diffuse incident radiation
    1978. 

  1978. 

  1979.              tauto=tausto+tauclf(i,j,k)
    1980. 

  1980.              ssato=(ssatau+ssacl(i,j,k)*tauclf(i,j,k))/tauto+1.0e-8
    1981. 

  1981.              ssato=min(ssato,0.999999)
    1982. 

  1982.              asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclf(i,j,k))/ &
    1983. 

  1983.                    (ssato*tauto)
    1984. 

  1984. 

  1985.              call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j))
    1986. 

  1986. 

  1987.             endif
    1988. 

  1988. 

  1989.            enddo
    1990. 

  1990.           enddo
    1991. 

  1991. 

  1992.            do j=1,n
    1993. 

  1993.             do i=1,m
    1994. 

  1994.                rr(i,j,k,1)=rr1t(i,j)
    1995. 

  1995.             enddo
    1996. 

  1996.            enddo
    1997. 

  1997.            do j=1,n
    1998. 

  1998.             do i=1,m
    1999. 

  1999.                tt(i,j,k,1)=tt1t(i,j)
    2000. 

  2000.             enddo
    2001. 

  2001.            enddo
    2002. 

  2002.            do j=1,n
    2003. 

  2003.             do i=1,m
    2004. 

  2004.                td(i,j,k,1)=td1t(i,j)
    2005. 

  2005.             enddo
    2006. 

  2006.            enddo
    2007. 

  2007.            do j=1,n
    2008. 

  2008.             do i=1,m
    2009. 

  2009.                rs(i,j,k,1)=rs1t(i,j)
    2010. 

  2010.             enddo
    2011. 

  2011.            enddo
    2012. 

  2012.            do j=1,n
    2013. 

  2013.             do i=1,m
    2014. 

  2014.                ts(i,j,k,1)=ts1t(i,j)
    2015. 

  2015.             enddo
    2016. 

  2016.            enddo
    2017. 

  2017.  
    2018. 

  2018.            do j=1,n
    2019. 

  2019.             do i=1,m
    2020. 

  2020.                rr(i,j,k,2)=rr2t(i,j)
    2021. 

  2021.             enddo
    2022. 

  2022.            enddo
    2023. 

  2023.            do j=1,n
    2024. 

  2024.             do i=1,m
    2025. 

  2025.                tt(i,j,k,2)=tt2t(i,j)
    2026. 

  2026.             enddo
    2027. 

  2027.            enddo
    2028. 

  2028.            do j=1,n
    2029. 

  2029.             do i=1,m
    2030. 

  2030.                td(i,j,k,2)=td2t(i,j)
    2031. 

  2031.             enddo
    2032. 

  2032.            enddo
    2033. 

  2033.            do j=1,n
    2034. 

  2034.             do i=1,m
    2035. 

  2035.                rs(i,j,k,2)=rs2t(i,j)
    2036. 

  2036.             enddo
    2037. 

  2037.            enddo
    2038. 

  2038.            do j=1,n
    2039. 

  2039.             do i=1,m
    2040. 

  2040.                ts(i,j,k,2)=ts2t(i,j)
    2041. 

  2041.             enddo
    2042. 

  2042.            enddo
    2043. 

  2043. 

  2044.  300  continue
    2045. 

  2045. 

  2046. !-----flux calculations
    2047. 

  2047. 

  2048.         call cldflx (m,n,np,ict,icb,overcast,cc,rr,tt,td,rs,ts, &
    2049. 

  2049.                      fclr,fall,fallu,falld,fsdir,fsdif)
    2050. 

  2050. 

  2051.        do k= 1, np+1
    2052. 

  2052.         do j= 1, n
    2053. 

  2053.          do i= 1, m
    2054. 

  2054.           flx(i,j,k) = flx(i,j,k)+fall(i,j,k)*hk(ib,ik)
    2055. 

  2055.           flxu(i,j,k) = flxu(i,j,k)+fallu(i,j,k)*hk(ib,ik)
    2056. 

  2056.           flxd(i,j,k) = flxd(i,j,k)+falld(i,j,k)*hk(ib,ik)
    2057. 

  2057.          enddo
    2058. 

  2058.         enddo
    2059. 

  2059.         do j= 1, n
    2060. 

  2060.          do i= 1, m
    2061. 

  2061.           flc(i,j,k) = flc(i,j,k)+fclr(i,j,k)*hk(ib,ik)
    2062. 

  2062.          enddo
    2063. 

  2063.         enddo
    2064. 

  2064.        enddo
    2065. 

  2065. 

  2066. !-----compute downward surface fluxes in the ir region
    2067. 

  2067. 

  2068.        do j= 1, n
    2069. 

  2069.         do i= 1, m
    2070. 

  2070.           fdirir(i,j) = fdirir(i,j)+fsdir(i,j)*hk(ib,ik)
    2071. 

  2071.           fdifir(i,j) = fdifir(i,j)+fsdif(i,j)*hk(ib,ik)
    2072. 

  2072.         enddo
    2073. 

  2073.        enddo
    2074. 

  2074. 

  2075.   200 continue
    2076. 

  2076.   100 continue
    2077. 

  2077.  
    2078. 

  2078.       end subroutine solir 
    2079. 

  2079. 

  2080. !********************************************************************
    2081. 

  2081. 

  2082.       subroutine cldscale (m,n,ndim,np,cosz,fcld,taucld,ict,icb,    &
    2083. 

  2083.                            cc,tauclb,tauclf)
    2084. 

  2084. 

  2085. !********************************************************************
    2086. 

  2086. !
    2087. 

  2087. !   This subroutine computes the high, middle, and
    2088. 

  2088. !    low cloud amounts and scales the cloud optical thickness.
    2089. 

  2089. !
    2090. 

  2090. !   To simplify calculations in a cloudy atmosphere, clouds are
    2091. 

  2091. !    grouped into high, middle and low clouds separated by the levels
    2092. 

  2092. !    ict and icb (level 1 is the top of the model atmosphere).
    2093. 

  2093. !
    2094. 

  2094. !   Within each of the three groups, clouds are assumed maximally
    2095. 

  2095. !    overlapped, and the cloud cover (cc) of a group is the maximum
    2096. 

  2096. !    cloud cover of all the layers in the group.  The optical thickness
    2097. 

  2097. !    (taucld) of a given layer is then scaled to new values (tauclb and
    2098. 

  2098. !    tauclf) so that the layer reflectance corresponding to the cloud
    2099. 

  2099. !    cover cc is the same as the original reflectance with optical
    2100. 

  2100. !    thickness taucld and cloud cover fcld.
    2101. 

  2101. !
    2102. 

  2102. !---input parameters
    2103. 

  2103. !
    2104. 

  2104. !    number of grid intervals in zonal direction (m)
    2105. 

  2105. !    number of grid intervals in meridional direction (n)
    2106. 

  2106. !    maximum number of grid intervals in meridional direction (ndim)
    2107. 

  2107. !    number of atmospheric layers (np)
    2108. 

  2108. !    cosine of the solar zenith angle (cosz)
    2109. 

  2109. !    fractional cloud cover (fcld)
    2110. 

  2110. !    cloud optical thickness (taucld)
    2111. 

  2111. !    index separating high and middle clouds (ict)
    2112. 

  2112. !    index separating middle and low clouds (icb)
    2113. 

  2113. !
    2114. 

  2114. !---output parameters
    2115. 

  2115. !
    2116. 

  2116. !    fractional cover of high, middle, and low clouds (cc)
    2117. 

  2117. !    scaled cloud optical thickness for beam radiation (tauclb)
    2118. 

  2118. !    scaled cloud optical thickness for diffuse radiation (tauclf)
    2119. 

  2119. !
    2120. 

  2120. !********************************************************************
    2121. 

  2121.       implicit none
    2122. 

  2122. !********************************************************************
    2123. 

  2123. 

  2124. !-----input parameters
    2125. 

  2125. 

  2126.       integer m,n,ndim,np
    2127. 

  2127.       integer ict(m,ndim),icb(m,ndim)
    2128. 

  2128.       real cosz(m,ndim),fcld(m,ndim,np),taucld(m,ndim,np,2)
    2129. 

  2129. 

  2130. !-----output parameters
    2131. 

  2131. 

  2132.       real cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)
    2133. 

  2133. 

  2134. !-----temporary variables
    2135. 

  2135. 

  2136.       integer i,j,k,im,it,ia,kk
    2137. 

  2137.       real  fm,ft,fa,xai,taux
    2138. 

  2138. 

  2139. !-----pre-computed table
    2140. 

  2140. 

  2141.       integer   nm,nt,na
    2142. 

  2142.       parameter (nm=11,nt=9,na=11) 
    2143. 

  2143.       real  dm,dt,da,t1,caib(nm,nt,na),caif(nt,na)
    2144. 

  2144.       parameter (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031)
    2145. 

  2145. 

  2146. !-----include the pre-computed table of mcai for scaling the cloud optical
    2147. 

  2147. !     thickness under the assumption that clouds are maximally overlapped
    2148. 

  2148. !
    2149. 

  2149. !     caib is for scaling the cloud optical thickness for direct radiation
    2150. 

  2150. !     caif is for scaling the cloud optical thickness for diffuse radiation
    2151. 

  2151. 

  2152. 

  2153.       data ((caib(1,i,j),j=1,11),i=1,9)/  &
    2154. 

  2154.        .000,0.068,0.140,0.216,0.298,0.385,0.481,0.586,0.705,0.840,1.000, &
    2155. 

  2155.        .000,0.052,0.106,0.166,0.230,0.302,0.383,0.478,0.595,0.752,1.000, &
    2156. 

  2156.        .000,0.038,0.078,0.120,0.166,0.218,0.276,0.346,0.438,0.582,1.000, &
    2157. 

  2157.        .000,0.030,0.060,0.092,0.126,0.164,0.206,0.255,0.322,0.442,1.000, &
    2158. 

  2158.        .000,0.025,0.051,0.078,0.106,0.136,0.170,0.209,0.266,0.462,1.000, &
    2159. 

  2159.        .000,0.023,0.046,0.070,0.095,0.122,0.150,0.187,0.278,0.577,1.000, &
    2160. 

  2160.        .000,0.022,0.043,0.066,0.089,0.114,0.141,0.187,0.354,0.603,1.000, &
    2161. 

  2161.        .000,0.021,0.042,0.063,0.086,0.108,0.135,0.214,0.349,0.565,1.000, &
    2162. 

  2162.        .000,0.021,0.041,0.062,0.083,0.105,0.134,0.202,0.302,0.479,1.000/
    2163. 

  2163.       data ((caib(2,i,j),j=1,11),i=1,9)/ &
    2164. 

  2164.        .000,0.088,0.179,0.272,0.367,0.465,0.566,0.669,0.776,0.886,1.000, &
    2165. 

  2165.        .000,0.079,0.161,0.247,0.337,0.431,0.531,0.637,0.749,0.870,1.000, &
    2166. 

  2166.        .000,0.065,0.134,0.207,0.286,0.372,0.466,0.572,0.692,0.831,1.000, &
    2167. 

  2167.        .000,0.049,0.102,0.158,0.221,0.290,0.370,0.465,0.583,0.745,1.000, &
    2168. 

  2168.        .000,0.037,0.076,0.118,0.165,0.217,0.278,0.354,0.459,0.638,1.000, &
    2169. 

  2169.        .000,0.030,0.061,0.094,0.130,0.171,0.221,0.286,0.398,0.631,1.000, &
    2170. 

  2170.        .000,0.026,0.052,0.081,0.111,0.146,0.189,0.259,0.407,0.643,1.000, &
    2171. 

  2171.        .000,0.023,0.047,0.072,0.098,0.129,0.170,0.250,0.387,0.598,1.000, &
    2172. 

  2172.        .000,0.022,0.044,0.066,0.090,0.118,0.156,0.224,0.328,0.508,1.000/
    2173. 

  2173.       data ((caib(3,i,j),j=1,11),i=1,9)/ &
    2174. 

  2174.        .000,0.094,0.189,0.285,0.383,0.482,0.582,0.685,0.788,0.894,1.000, &
    2175. 

  2175.        .000,0.088,0.178,0.271,0.366,0.465,0.565,0.669,0.776,0.886,1.000, &
    2176. 

  2176.        .000,0.079,0.161,0.247,0.337,0.431,0.531,0.637,0.750,0.870,1.000, &
    2177. 

  2177.        .000,0.066,0.134,0.209,0.289,0.375,0.470,0.577,0.697,0.835,1.000, &
    2178. 

  2178.        .000,0.050,0.104,0.163,0.227,0.300,0.383,0.483,0.606,0.770,1.000, &
    2179. 

  2179.        .000,0.038,0.080,0.125,0.175,0.233,0.302,0.391,0.518,0.710,1.000, &
    2180. 

  2180.        .000,0.031,0.064,0.100,0.141,0.188,0.249,0.336,0.476,0.689,1.000, &
    2181. 

  2181.        .000,0.026,0.054,0.084,0.118,0.158,0.213,0.298,0.433,0.638,1.000, &
    2182. 

  2182.        .000,0.023,0.048,0.074,0.102,0.136,0.182,0.254,0.360,0.542,1.000/
    2183. 

  2183.       data ((caib(4,i,j),j=1,11),i=1,9)/ &
    2184. 

  2184.        .000,0.096,0.193,0.290,0.389,0.488,0.589,0.690,0.792,0.896,1.000, &
    2185. 

  2185.        .000,0.092,0.186,0.281,0.378,0.477,0.578,0.680,0.785,0.891,1.000, &
    2186. 

  2186.        .000,0.086,0.174,0.264,0.358,0.455,0.556,0.660,0.769,0.882,1.000, &
    2187. 

  2187.        .000,0.074,0.153,0.235,0.323,0.416,0.514,0.622,0.737,0.862,1.000, &
    2188. 

  2188.        .000,0.061,0.126,0.195,0.271,0.355,0.449,0.555,0.678,0.823,1.000, &
    2189. 

  2189.        .000,0.047,0.098,0.153,0.215,0.286,0.370,0.471,0.600,0.770,1.000, &
    2190. 

  2190.        .000,0.037,0.077,0.120,0.170,0.230,0.303,0.401,0.537,0.729,1.000, &
    2191. 

  2191.        .000,0.030,0.062,0.098,0.138,0.187,0.252,0.343,0.476,0.673,1.000, &
    2192. 

  2192.        .000,0.026,0.053,0.082,0.114,0.154,0.207,0.282,0.391,0.574,1.000/
    2193. 

  2193.       data ((caib(5,i,j),j=1,11),i=1,9)/ &
    2194. 

  2194.        .000,0.097,0.194,0.293,0.392,0.492,0.592,0.693,0.794,0.897,1.000, &
    2195. 

  2195.        .000,0.094,0.190,0.286,0.384,0.483,0.584,0.686,0.789,0.894,1.000, &
    2196. 

  2196.        .000,0.090,0.181,0.274,0.370,0.468,0.569,0.672,0.778,0.887,1.000, &
    2197. 

  2197.        .000,0.081,0.165,0.252,0.343,0.439,0.539,0.645,0.757,0.874,1.000, &
    2198. 

  2198.        .000,0.069,0.142,0.218,0.302,0.392,0.490,0.598,0.717,0.850,1.000, &
    2199. 

  2199.        .000,0.054,0.114,0.178,0.250,0.330,0.422,0.529,0.656,0.810,1.000, &
    2200. 

  2200.        .000,0.042,0.090,0.141,0.200,0.269,0.351,0.455,0.589,0.764,1.000, &
    2201. 

  2201.        .000,0.034,0.070,0.112,0.159,0.217,0.289,0.384,0.515,0.703,1.000, &
    2202. 

  2202.        .000,0.028,0.058,0.090,0.128,0.174,0.231,0.309,0.420,0.602,1.000/
    2203. 

  2203.       data ((caib(6,i,j),j=1,11),i=1,9)/ &
    2204. 

  2204.        .000,0.098,0.196,0.295,0.394,0.494,0.594,0.695,0.796,0.898,1.000, &
    2205. 

  2205.        .000,0.096,0.193,0.290,0.389,0.488,0.588,0.690,0.792,0.895,1.000, &
    2206. 

  2206.        .000,0.092,0.186,0.281,0.378,0.477,0.577,0.680,0.784,0.891,1.000, &
    2207. 

  2207.        .000,0.086,0.174,0.264,0.358,0.455,0.556,0.661,0.769,0.882,1.000, &
    2208. 

  2208.        .000,0.075,0.154,0.237,0.325,0.419,0.518,0.626,0.741,0.865,1.000, &
    2209. 

  2209.        .000,0.062,0.129,0.201,0.279,0.366,0.462,0.571,0.694,0.836,1.000, &
    2210. 

  2210.        .000,0.049,0.102,0.162,0.229,0.305,0.394,0.501,0.631,0.793,1.000, &
    2211. 

  2211.        .000,0.038,0.080,0.127,0.182,0.245,0.323,0.422,0.550,0.730,1.000, &
    2212. 

  2212.        .000,0.030,0.064,0.100,0.142,0.192,0.254,0.334,0.448,0.627,1.000/
    2213. 

  2213.       data ((caib(7,i,j),j=1,11),i=1,9)/ &
    2214. 

  2214.        .000,0.098,0.198,0.296,0.396,0.496,0.596,0.696,0.797,0.898,1.000, &
    2215. 

  2215.        .000,0.097,0.194,0.293,0.392,0.491,0.591,0.693,0.794,0.897,1.000, &
    2216. 

  2216.        .000,0.094,0.190,0.286,0.384,0.483,0.583,0.686,0.789,0.894,1.000, &
    2217. 

  2217.        .000,0.089,0.180,0.274,0.369,0.467,0.568,0.672,0.778,0.887,1.000, &
    2218. 

  2218.        .000,0.081,0.165,0.252,0.344,0.440,0.541,0.646,0.758,0.875,1.000, &
    2219. 

  2219.        .000,0.069,0.142,0.221,0.306,0.397,0.496,0.604,0.722,0.854,1.000, &
    2220. 

  2220.        .000,0.056,0.116,0.182,0.256,0.338,0.432,0.540,0.666,0.816,1.000, &
    2221. 

  2221.        .000,0.043,0.090,0.143,0.203,0.273,0.355,0.455,0.583,0.754,1.000, &
    2222. 

  2222.        .000,0.034,0.070,0.111,0.157,0.210,0.276,0.359,0.474,0.650,1.000/
    2223. 

  2223.       data ((caib(8,i,j),j=1,11),i=1,9)/ &
    2224. 

  2224.        .000,0.099,0.198,0.298,0.398,0.497,0.598,0.698,0.798,0.899,1.000, &
    2225. 

  2225.        .000,0.098,0.196,0.295,0.394,0.494,0.594,0.695,0.796,0.898,1.000, &
    2226. 

  2226.        .000,0.096,0.193,0.290,0.390,0.489,0.589,0.690,0.793,0.896,1.000, &
    2227. 

  2227.        .000,0.093,0.186,0.282,0.379,0.478,0.578,0.681,0.786,0.892,1.000, &
    2228. 

  2228.        .000,0.086,0.175,0.266,0.361,0.458,0.558,0.663,0.771,0.883,1.000, &
    2229. 

  2229.        .000,0.076,0.156,0.240,0.330,0.423,0.523,0.630,0.744,0.867,1.000, &
    2230. 

  2230.        .000,0.063,0.130,0.203,0.282,0.369,0.465,0.572,0.694,0.834,1.000, &
    2231. 

  2231.        .000,0.049,0.102,0.161,0.226,0.299,0.385,0.486,0.611,0.774,1.000, &
    2232. 

  2232.        .000,0.038,0.078,0.122,0.172,0.229,0.297,0.382,0.498,0.672,1.000/
    2233. 

  2233.       data ((caib(9,i,j),j=1,11),i=1,9)/ &
    2234. 

  2234.        .000,0.099,0.199,0.298,0.398,0.498,0.598,0.699,0.799,0.899,1.000, &
    2235. 

  2235.        .000,0.099,0.198,0.298,0.398,0.497,0.598,0.698,0.798,0.899,1.000, &
    2236. 

  2236.        .000,0.098,0.196,0.295,0.394,0.494,0.594,0.695,0.796,0.898,1.000, &
    2237. 

  2237.        .000,0.096,0.193,0.290,0.389,0.488,0.588,0.690,0.792,0.895,1.000, &
    2238. 

  2238.        .000,0.092,0.185,0.280,0.376,0.474,0.575,0.678,0.782,0.890,1.000, &
    2239. 

  2239.        .000,0.084,0.170,0.259,0.351,0.447,0.547,0.652,0.762,0.878,1.000, &
    2240. 

  2240.        .000,0.071,0.146,0.224,0.308,0.398,0.494,0.601,0.718,0.850,1.000, &
    2241. 

  2241.        .000,0.056,0.114,0.178,0.248,0.325,0.412,0.514,0.638,0.793,1.000, &
    2242. 

  2242.        .000,0.042,0.086,0.134,0.186,0.246,0.318,0.405,0.521,0.691,1.000/
    2243. 

  2243.       data ((caib(10,i,j),j=1,11),i=1,9)/ &
    2244. 

  2244.        .000,0.100,0.200,0.300,0.400,0.500,0.600,0.700,0.800,0.900,1.000, &
    2245. 

  2245.        .000,0.100,0.200,0.300,0.400,0.500,0.600,0.700,0.800,0.900,1.000, &
    2246. 

  2246.        .000,0.100,0.200,0.300,0.400,0.500,0.600,0.700,0.800,0.900,1.000, &
    2247. 

  2247.        .000,0.100,0.199,0.298,0.398,0.498,0.598,0.698,0.798,0.899,1.000, &
    2248. 

  2248.        .000,0.098,0.196,0.294,0.392,0.491,0.590,0.691,0.793,0.896,1.000, &
    2249. 

  2249.        .000,0.092,0.185,0.278,0.374,0.470,0.570,0.671,0.777,0.886,1.000, &
    2250. 

  2250.        .000,0.081,0.162,0.246,0.333,0.424,0.521,0.625,0.738,0.862,1.000, &
    2251. 

  2251.        .000,0.063,0.128,0.196,0.270,0.349,0.438,0.540,0.661,0.809,1.000, &
    2252. 

  2252.        .000,0.046,0.094,0.146,0.202,0.264,0.337,0.426,0.542,0.710,1.000/ 
    2253. 

  2253.       data ((caib(11,i,j),j=1,11),i=1,9)/ &
    2254. 

  2254.        .000,0.101,0.202,0.302,0.402,0.502,0.602,0.702,0.802,0.901,1.000, &
    2255. 

  2255.        .000,0.102,0.202,0.303,0.404,0.504,0.604,0.703,0.802,0.902,1.000, &
    2256. 

  2256.        .000,0.102,0.205,0.306,0.406,0.506,0.606,0.706,0.804,0.902,1.000, &
    2257. 

  2257.        .000,0.104,0.207,0.309,0.410,0.510,0.609,0.707,0.805,0.902,1.000, &
    2258. 

  2258.        .000,0.106,0.208,0.309,0.409,0.508,0.606,0.705,0.803,0.902,1.000, &
    2259. 

  2259.        .000,0.102,0.202,0.298,0.395,0.493,0.590,0.690,0.790,0.894,1.000, &
    2260. 

  2260.        .000,0.091,0.179,0.267,0.357,0.449,0.545,0.647,0.755,0.872,1.000, &
    2261. 

  2261.        .000,0.073,0.142,0.214,0.290,0.372,0.462,0.563,0.681,0.822,1.000, &
    2262. 

  2262.        .000,0.053,0.104,0.158,0.217,0.281,0.356,0.446,0.562,0.726,1.000/
    2263. 

  2263.       data ((caif(i,j),j=1,11),i=1,9)/ &
    2264. 

  2264.        .000,0.099,0.198,0.297,0.397,0.496,0.597,0.697,0.798,0.899,1.000, &
    2265. 

  2265.        .000,0.098,0.196,0.294,0.394,0.494,0.594,0.694,0.796,0.898,1.000, &
    2266. 

  2266.        .000,0.096,0.192,0.290,0.388,0.487,0.587,0.689,0.792,0.895,1.000, &
    2267. 

  2267.        .000,0.092,0.185,0.280,0.376,0.476,0.576,0.678,0.783,0.890,1.000, &
    2268. 

  2268.        .000,0.085,0.173,0.263,0.357,0.454,0.555,0.659,0.768,0.881,1.000, &
    2269. 

  2269.        .000,0.076,0.154,0.237,0.324,0.418,0.517,0.624,0.738,0.864,1.000, &
    2270. 

  2270.        .000,0.063,0.131,0.203,0.281,0.366,0.461,0.567,0.688,0.830,1.000, &
    2271. 

  2271.        .000,0.052,0.107,0.166,0.232,0.305,0.389,0.488,0.610,0.770,1.000, &
    2272. 

  2272.        .000,0.043,0.088,0.136,0.189,0.248,0.317,0.400,0.510,0.675,1.000/
    2273. 

  2273. 

  2274. !-----clouds within each of the high, middle, and low clouds are assumed
    2275. 

  2275. !     to be maximally overlapped, and the cloud cover (cc) for a group
    2276. 

  2276. !     (high, middle, or low) is the maximum cloud cover of all the layers
    2277. 

  2277. !     within a group
    2278. 

  2278. 

  2279.       do j=1,n
    2280. 

  2280.        do i=1,m
    2281. 

  2281.          cc(i,j,1)=0.0
    2282. 

  2282.          cc(i,j,2)=0.0
    2283. 

  2283.          cc(i,j,3)=0.0
    2284. 

  2284.        enddo
    2285. 

  2285.       enddo
    2286. 

  2286.       do j=1,n
    2287. 

  2287.        do i=1,m
    2288. 

  2288.         do k=1,ict(i,j)-1
    2289. 

  2289.           cc(i,j,1)=max(cc(i,j,1),fcld(i,j,k))
    2290. 

  2290.         enddo
    2291. 

  2291.        enddo
    2292. 

  2292.       enddo
    2293. 

  2293. 

  2294.       do j=1,n
    2295. 

  2295.        do i=1,m
    2296. 

  2296.         do k=ict(i,j),icb(i,j)-1
    2297. 

  2297.           cc(i,j,2)=max(cc(i,j,2),fcld(i,j,k))
    2298. 

  2298.         enddo
    2299. 

  2299.        enddo
    2300. 

  2300.       enddo
    2301. 

  2301. 

  2302.       do j=1,n
    2303. 

  2303.        do i=1,m
    2304. 

  2304.         do k=icb(i,j),np
    2305. 

  2305.           cc(i,j,3)=max(cc(i,j,3),fcld(i,j,k))
    2306. 

  2306.         enddo
    2307. 

  2307.        enddo
    2308. 

  2308.       enddo
    2309. 

  2309. 

  2310. !-----scale the cloud optical thickness.
    2311. 

  2311. !     taucld(i,j,k,1) is the optical thickness for ice particles, and
    2312. 

  2312. !     taucld(i,j,k,2) is the optical thickness for liquid particles.
    2313. 

  2313.       
    2314. 

  2314.       do j=1,n
    2315. 

  2315.        do i=1,m
    2316. 

  2316. 

  2317.         do k=1,np
    2318. 

  2318. 

  2319.          if(k.lt.ict(i,j)) then
    2320. 

  2320.             kk=1
    2321. 

  2321.          elseif(k.ge.ict(i,j) .and. k.lt.icb(i,j)) then
    2322. 

  2322.             kk=2
    2323. 

  2323.          else
    2324. 

  2324.             kk=3
    2325. 

  2325.          endif
    2326. 

  2326. 

  2327.          tauclb(i,j,k) = 0.0
    2328. 

  2328.          tauclf(i,j,k) = 0.0
    2329. 

  2329. 

  2330.          taux=taucld(i,j,k,1)+taucld(i,j,k,2)
    2331. 

  2331.          if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then
    2332. 

  2332. 

  2333. !-----normalize cloud cover
    2334. 

  2334. 

  2335.            fa=fcld(i,j,k)/cc(i,j,kk)
    2336. 

  2336. 

  2337. !-----table look-up
    2338. 

  2338. 

  2339.            taux=min(taux,32.)
    2340. 

  2340. 

  2341.            fm=cosz(i,j)/dm
    2342. 

  2342.            ft=(log10(taux)-t1)/dt
    2343. 

  2343.            fa=fa/da
    2344. 

  2344.  
    2345. 

  2345.            im=int(fm+1.5)
    2346. 

  2346.            it=int(ft+1.5)
    2347. 

  2347.            ia=int(fa+1.5)
    2348. 

  2348.   
    2349. 

  2349.            im=max(im,2)
    2350. 

  2350.            it=max(it,2)
    2351. 

  2351.            ia=max(ia,2)
    2352. 

  2352.      
    2353. 

  2353.            im=min(im,nm-1)
    2354. 

  2354.            it=min(it,nt-1)
    2355. 

  2355.            ia=min(ia,na-1)
    2356. 

  2356. 

  2357.            fm=fm-float(im-1)
    2358. 

  2358.            ft=ft-float(it-1)
    2359. 

  2359.            fa=fa-float(ia-1)
    2360. 

  2360. 

  2361. !-----scale cloud optical thickness for beam radiation.
    2362. 

  2362. !     the scaling factor, xai, is a function of the solar zenith
    2363. 

  2363. !     angle, optical thickness, and cloud cover.
    2364. 

  2364.  
    2365. 

  2365.            xai=    (-caib(im-1,it,ia)*(1.-fm)+ &
    2366. 

  2366.             caib(im+1,it,ia)*(1.+fm))*fm*.5+caib(im,it,ia)*(1.-fm*fm)
    2367. 

  2367.          
    2368. 

  2368.            xai=xai+(-caib(im,it-1,ia)*(1.-ft)+ &
    2369. 

  2369.             caib(im,it+1,ia)*(1.+ft))*ft*.5+caib(im,it,ia)*(1.-ft*ft)
    2370. 

  2370. 

  2371.            xai=xai+(-caib(im,it,ia-1)*(1.-fa)+ &
    2372. 

  2372.            caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa)
    2373. 

  2373. 

  2374.            xai= xai-2.*caib(im,it,ia)
    2375. 

  2375.            xai=max(xai,0.0)
    2376. 

  2376.      
    2377. 

  2377.            tauclb(i,j,k) = taux*xai
    2378. 

  2378. 

  2379. !-----scale cloud optical thickness for diffuse radiation.
    2380. 

  2380. !     the scaling factor, xai, is a function of the cloud optical
    2381. 

  2381. !     thickness and cover but not the solar zenith angle.
    2382. 

  2382. 

  2383.            xai=    (-caif(it-1,ia)*(1.-ft)+ &
    2384. 

  2384.             caif(it+1,ia)*(1.+ft))*ft*.5+caif(it,ia)*(1.-ft*ft)
    2385. 

  2385. 

  2386.            xai=xai+(-caif(it,ia-1)*(1.-fa)+ &
    2387. 

  2387.             caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa)
    2388. 

  2388. 

  2389.            xai= xai-caif(it,ia)
    2390. 

  2390.            xai=max(xai,0.0)
    2391. 

  2391.      
    2392. 

  2392.            tauclf(i,j,k) = taux*xai
    2393. 

  2393. 

  2394.          endif
    2395. 

  2395. 

  2396.         enddo
    2397. 

  2397.        enddo
    2398. 

  2398.       enddo
    2399. 

  2399. 

  2400.       end subroutine cldscale
    2401. 

  2401. 

  2402. !*********************************************************************
    2403. 

  2403. 

  2404.       subroutine deledd(tau,ssc,g0,csm,rr,tt,td)
    2405. 

  2405. 

  2406. !*********************************************************************
    2407. 

  2407. !
    2408. 

  2408. !-----uses the delta-eddington approximation to compute the
    2409. 

  2409. !     bulk scattering properties of a single layer
    2410. 

  2410. !     coded following King and Harshvardhan (JAS, 1986)
    2411. 

  2411. !
    2412. 

  2412. !  inputs:
    2413. 

  2413. !
    2414. 

  2414. !     tau: the effective optical thickness
    2415. 

  2415. !     ssc: the effective single scattering albedo
    2416. 

  2416. !     g0:  the effective asymmetry factor
    2417. 

  2417. !     csm: the effective secant of the zenith angle
    2418. 

  2418. !
    2419. 

  2419. !  outputs:
    2420. 

  2420. !
    2421. 

  2421. !     rr: the layer reflection of the direct beam
    2422. 

  2422. !     tt: the layer diffuse transmission of the direct beam
    2423. 

  2423. !     td: the layer direct transmission of the direct beam
    2424. 

  2424. !
    2425. 

  2425. !*********************************************************************
    2426. 

  2426.       implicit none
    2427. 

  2427. !*********************************************************************
    2428. 

  2428. 

  2429.       real zero,one,two,three,four,fourth,seven,thresh
    2430. 

  2430.       parameter (one =1., three=3.)
    2431. 

  2431.       parameter (two =2., seven=7.)
    2432. 

  2432.       parameter (four=4., fourth=.25)
    2433. 

  2433.       parameter (zero=0., thresh=1.e-8)
    2434. 

  2434. 

  2435. !-----input parameters
    2436. 

  2436.       real tau,ssc,g0,csm
    2437. 

  2437. 

  2438. !-----output parameters
    2439. 

  2439.       real rr,tt,td
    2440. 

  2440. 

  2441. !-----temporary parameters
    2442. 

  2442. 

  2443.       real zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2, &
    2444. 

  2444.            all,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4
    2445. 

  2445.  
    2446. 

  2446. !---------------------------------------------------------------------
    2447. 

  2447. 

  2448.                 zth = one / csm
    2449. 

  2449.  
    2450. 

  2450. !  delta-eddington scaling of single scattering albedo,
    2451. 

  2451. !  optical thickness, and asymmetry factor,
    2452. 

  2452. !  K & H eqs(27-29)
    2453. 

  2453. 

  2454.                 ff  = g0*g0
    2455. 

  2455.                 xx  = one-ff*ssc
    2456. 

  2456.                 taup= tau*xx
    2457. 

  2457.                 sscp= ssc*(one-ff)/xx
    2458. 

  2458.                 gp  = g0/(one+g0)
    2459. 

  2459.  
    2460. 

  2460. !  gamma1, gamma2, and gamma3. see table 2 and eq(26) K & H
    2461. 

  2461. !  ssc and gp are the d-s single scattering
    2462. 

  2462. !  albedo and asymmetry factor.
    2463. 

  2463. 

  2464.                 xx  =  three*gp 
    2465. 

  2465.                 gm1 =  (seven - sscp*(four+xx))*fourth
    2466. 

  2466.                 gm2 = -(one   - sscp*(four-xx))*fourth
    2467. 

  2467.  
    2468. 

  2468. !  akk is k as defined in eq(25) of K & H
    2469. 

  2469.  
    2470. 

  2470.                 akk = sqrt((gm1+gm2)*(gm1-gm2))
    2471. 

  2471.  
    2472. 

  2472.                 xx  = akk * zth
    2473. 

  2473.                 st7 = one - xx
    2474. 

  2474.                 st8 = one + xx
    2475. 

  2475.                 st3 = st7 * st8
    2476. 

  2476. 

  2477.                 if (abs(st3) .lt. thresh) then
    2478. 

  2478.                     zth = zth + 0.001
    2479. 

  2479.                     xx  = akk * zth
    2480. 

  2480.                     st7 = one - xx
    2481. 

  2481.                     st8 = one + xx
    2482. 

  2482.                     st3 = st7 * st8
    2483. 

  2483.                 endif
    2484. 

  2484. 

  2485. !  extinction of the direct beam transmission
    2486. 

  2486.  
    2487. 

  2487.                 td  = exp(-taup/zth)
    2488. 

  2488. 

  2489. !  alf1 and alf2 are alpha1 and alpha2 from eqs (23) & (24) of K & H
    2490. 

  2490.  
    2491. 

  2491.                 gm3  = (two - zth*three*gp)*fourth
    2492. 

  2492.                 xx   = gm1 - gm2
    2493. 

  2493.                 alf1 = gm1 - gm3 * xx
    2494. 

  2494.                 alf2 = gm2 + gm3 * xx
    2495. 

  2495.  
    2496. 

  2496. !  all is last term in eq(21) of K & H
    2497. 

  2497. !  bll is last term in eq(22) of K & H
    2498. 

  2498.  
    2499. 

  2499.                 xx  = akk * two
    2500. 

  2500.                 all = (gm3 - alf2 * zth    )*xx*td
    2501. 

  2501.                 bll = (one - gm3 + alf1*zth)*xx
    2502. 

  2502.  
    2503. 

  2503.                 xx  = akk * gm3
    2504. 

  2504.                 cll = (alf2 + xx) * st7
    2505. 

  2505.                 dll = (alf2 - xx) * st8
    2506. 

  2506.  
    2507. 

  2507.                 xx  = akk * (one-gm3)
    2508. 

  2508.                 fll = (alf1 + xx) * st8
    2509. 

  2509.                 ell = (alf1 - xx) * st7
    2510. 

  2510.   
    2511. 

  2511.                 st2 = exp(-akk*taup)
    2512. 

  2512.                 st4 = st2 * st2
    2513. 

  2513.  
    2514. 

  2514.                 st1 =  sscp / ((akk+gm1 + (akk-gm1)*st4) * st3)
    2515. 

  2515.  
    2516. 

  2516. !  rr is r-hat of eq(21) of K & H
    2517. 

  2517. !  tt is diffuse part of t-hat of eq(22) of K & H
    2518. 

  2518.  
    2519. 

  2519.                 rr =   ( cll-dll*st4    -all*st2)*st1
    2520. 

  2520.                 tt = - ((fll-ell*st4)*td-bll*st2)*st1
    2521. 

  2521.  
    2522. 

  2522.                 rr = max(rr,zero)
    2523. 

  2523.                 tt = max(tt,zero)
    2524. 

  2524. 

  2525.       end subroutine deledd
    2526. 

  2526. 

  2527. !*********************************************************************
    2528. 

  2528. 

  2529.       subroutine sagpol(tau,ssc,g0,rll,tll)
    2530. 

  2530. 

  2531. !*********************************************************************
    2532. 

  2532. !-----transmittance (tll) and reflectance (rll) of diffuse radiation
    2533. 

  2533. !     follows Sagan and Pollock (JGR, 1967).
    2534. 

  2534. !     also, eq.(31) of Lacis and Hansen (JAS, 1974).
    2535. 

  2535. !
    2536. 

  2536. !-----input parameters:
    2537. 

  2537. !
    2538. 

  2538. !      tau: the effective optical thickness
    2539. 

  2539. !      ssc: the effective single scattering albedo
    2540. 

  2540. !      g0:  the effective asymmetry factor
    2541. 

  2541. !
    2542. 

  2542. !-----output parameters:
    2543. 

  2543. !
    2544. 

  2544. !      rll: the layer reflection of diffuse radiation
    2545. 

  2545. !      tll: the layer transmission of diffuse radiation
    2546. 

  2546. !
    2547. 

  2547. !*********************************************************************
    2548. 

  2548.       implicit none
    2549. 

  2549. !*********************************************************************
    2550. 

  2550. 

  2551.       real one,three,four
    2552. 

  2552.       parameter (one=1., three=3., four=4.)
    2553. 

  2553. 

  2554. !-----output parameters:
    2555. 

  2555. 

  2556.       real tau,ssc,g0
    2557. 

  2557. 

  2558. !-----output parameters:
    2559. 

  2559. 

  2560.       real rll,tll
    2561. 

  2561. 

  2562. !-----temporary arrays
    2563. 

  2563. 

  2564.       real xx,uuu,ttt,emt,up1,um1,st1
    2565. 

  2565. 

  2566.              xx  = one-ssc*g0
    2567. 

  2567.              uuu = sqrt( xx/(one-ssc))
    2568. 

  2568.              ttt = sqrt( xx*(one-ssc)*three )*tau
    2569. 

  2569.              emt = exp(-ttt)
    2570. 

  2570.              up1 = uuu + one
    2571. 

  2571.              um1 = uuu - one
    2572. 

  2572.              xx  = um1*emt
    2573. 

  2573.              st1 = one / ((up1+xx) * (up1-xx))
    2574. 

  2574.              rll = up1*um1*(one-emt*emt)*st1
    2575. 

  2575.              tll = uuu*four*emt         *st1
    2576. 

  2576. 

  2577.       end subroutine sagpol
    2578. 

  2578. 

  2579. !*******************************************************************
    2580. 

  2580. 

  2581.       subroutine cldflx (m,n,np,ict,icb,overcast,cc,rr,tt,td,rs,ts,&
    2582. 

  2582.                  fclr,fall,fallu,falld,fsdir,fsdif)
    2583. 

  2583. 

  2584. !*******************************************************************
    2585. 

  2585. !  compute upward and downward fluxes using a two-stream adding method
    2586. 

  2586. !  following equations (3)-(5) of Chou (1992, JAS).
    2587. 

  2587. !
    2588. 

  2588. !  clouds are grouped into high, middle, and low clouds which are 
    2589. 

  2589. !  assumed randomly overlapped. It involves eight sets of calculations.
    2590. 

  2590. !  In each set of calculations, each atmospheric layer is homogeneous,
    2591. 

  2591. !  either totally filled with clouds or without clouds.
    2592. 

  2592. 

  2593. !  input parameters:
    2594. 

  2594. !
    2595. 

  2595. !     m:   number of soundings in zonal direction
    2596. 

  2596. !     n:   number of soundings in meridional direction
    2597. 

  2597. !     np:  number of atmospheric layers
    2598. 

  2598. !     ict: the level separating high and middle clouds
    2599. 

  2599. !     icb: the level separating middle and low clouds
    2600. 

  2600. !     cc:  effective cloud covers for high, middle and low clouds
    2601. 

  2601. !     tt:  diffuse transmission of a layer illuminated by beam radiation
    2602. 

  2602. !     td:  direct beam tranmssion
    2603. 

  2603. !     ts:  transmission of a layer illuminated by diffuse radiation
    2604. 

  2604. !     rr:  reflection of a layer illuminated by beam radiation
    2605. 

  2605. !     rs:  reflection of a layer illuminated by diffuse radiation
    2606. 

  2606. !
    2607. 

  2607. !  output parameters:
    2608. 

  2608. !
    2609. 

  2609. !     fclr:  clear-sky flux (downward minus upward)
    2610. 

  2610. !     fall:  all-sky flux (downward minus upward)
    2611. 

  2611. !     fsdir: surface direct downward flux
    2612. 

  2612. !     fsdif: surface diffuse downward flux
    2613. 

  2613. !
    2614. 

  2614. !*********************************************************************c
    2615. 

  2615.       implicit none
    2616. 

  2616. !*********************************************************************c
    2617. 

  2617. 

  2618. !-----input parameters
    2619. 

  2619. 

  2620.       integer m,n,np
    2621. 

  2621.       integer ict(m,n),icb(m,n)
    2622. 

  2622. 

  2623.       real rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2)
    2624. 

  2624.       real rs(m,n,np+1,2),ts(m,n,np+1,2)
    2625. 

  2625.       real cc(m,n,3)
    2626. 

  2626.       logical overcast
    2627. 

  2627. 

  2628. !-----temporary array
    2629. 

  2629. 

  2630.       integer i,j,k,ih,im,is,itm
    2631. 

  2631.       real rra(m,n,np+1,2,2),tta(m,n,np+1,2,2),tda(m,n,np+1,2,2)
    2632. 

  2632.       real rsa(m,n,np+1,2,2),rxa(m,n,np+1,2,2)
    2633. 

  2633.       real ch(m,n),cm(m,n),ct(m,n),flxdn(m,n,np+1)
    2634. 

  2634.       real flxdnu(m,n,np+1),flxdnd(m,n,np+1)
    2635. 

  2635.       real fdndir(m,n),fdndif(m,n),fupdif
    2636. 

  2636.       real denm,xx
    2637. 

  2637. 

  2638. !-----output parameters
    2639. 

  2639. 

  2640.       real fclr(m,n,np+1),fall(m,n,np+1)
    2641. 

  2641.       real fallu(m,n,np+1),falld(m,n,np+1)
    2642. 

  2642.       real fsdir(m,n),fsdif(m,n)
    2643. 

  2643. 

  2644. !-----initialize all-sky flux (fall) and surface downward fluxes
    2645. 

  2645. 

  2646.       do k=1,np+1
    2647. 

  2647.        do j=1,n
    2648. 

  2648.         do i=1,m
    2649. 

  2649.            fclr(i,j,k)=0.0
    2650. 

  2650.            fall(i,j,k)=0.0
    2651. 

  2651.            fallu(i,j,k)=0.0
    2652. 

  2652.            falld(i,j,k)=0.0
    2653. 

  2653.         enddo
    2654. 

  2654.        enddo
    2655. 

  2655.       enddo
    2656. 

  2656. 

  2657.        do j=1,n
    2658. 

  2658.         do i=1,m
    2659. 

  2659.            fsdir(i,j)=0.0
    2660. 

  2660.            fsdif(i,j)=0.0
    2661. 

  2661.         enddo
    2662. 

  2662.        enddo
    2663. 

  2663. 

  2664. !-----compute transmittances and reflectances for a composite of
    2665. 

  2665. !     layers. layers are added one at a time, going down from the top.
    2666. 

  2666. !     tda is the composite transmittance illuminated by beam radiation
    2667. 

  2667. !     tta is the composite diffuse transmittance illuminated by
    2668. 

  2668. !         beam radiation
    2669. 

  2669. !     rsa is the composite reflectance illuminated from below
    2670. 

  2670. !         by diffuse radiation
    2671. 

  2671. !     tta and rsa are computed from eqs. (4b) and (3b) of Chou
    2672. 

  2672. 

  2673.       itm=1
    2674. 

  2674. 

  2675. !-----if overcas.=.true., set itm=2, and only one set of fluxes is computed
    2676. 

  2676. 

  2677.       if (overcast) itm=2
    2678. 

  2678. 

  2679. !-----for high clouds. indices 1 and 2 denote clear and cloudy
    2680. 

  2680. !     situations, respectively.
    2681. 

  2681. 

  2682.       do 10 ih=itm,2
    2683. 

  2683. 

  2684.        do j= 1, n
    2685. 

  2685.         do i= 1, m
    2686. 

  2686.           tda(i,j,1,ih,1)=td(i,j,1,ih)
    2687. 

  2687.           tta(i,j,1,ih,1)=tt(i,j,1,ih)
    2688. 

  2688.           rsa(i,j,1,ih,1)=rs(i,j,1,ih)
    2689. 

  2689.           tda(i,j,1,ih,2)=td(i,j,1,ih)
    2690. 

  2690.           tta(i,j,1,ih,2)=tt(i,j,1,ih)
    2691. 

  2691.           rsa(i,j,1,ih,2)=rs(i,j,1,ih)
    2692. 

  2692.         enddo
    2693. 

  2693.        enddo
    2694. 

  2694. 

  2695.        do j= 1, n
    2696. 

  2696.         do i= 1, m
    2697. 

  2697.          do k= 2, ict(i,j)-1
    2698. 

  2698.           denm = ts(i,j,k,ih)/( 1.-rsa(i,j,k-1,ih,1)*rs(i,j,k,ih))
    2699. 

  2699.           tda(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*td(i,j,k,ih)
    2700. 

  2700.           tta(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*tt(i,j,k,ih) &
    2701. 

  2701.                         +(tda(i,j,k-1,ih,1)*rr(i,j,k,ih)  &
    2702. 

  2702.                         *rsa(i,j,k-1,ih,1)+tta(i,j,k-1,ih,1))*denm
    2703. 

  2703.           rsa(i,j,k,ih,1)= rs(i,j,k,ih)+ts(i,j,k,ih) &
    2704. 

  2704.                         *rsa(i,j,k-1,ih,1)*denm
    2705. 

  2705.           tda(i,j,k,ih,2)= tda(i,j,k,ih,1)
    2706. 

  2706.           tta(i,j,k,ih,2)= tta(i,j,k,ih,1)
    2707. 

  2707.           rsa(i,j,k,ih,2)= rsa(i,j,k,ih,1)
    2708. 

  2708.          enddo
    2709. 

  2709.         enddo
    2710. 

  2710.        enddo
    2711. 

  2711. 

  2712. !-----for middle clouds
    2713. 

  2713. 

  2714.       do 10 im=itm,2
    2715. 

  2715. 

  2716.       do j= 1, n
    2717. 

  2717.        do i= 1, m
    2718. 

  2718.         do k= ict(i,j), icb(i,j)-1
    2719. 

  2719.           denm = ts(i,j,k,im)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,im))
    2720. 

  2720.           tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,im)
    2721. 

  2721.           tta(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*tt(i,j,k,im) &
    2722. 

  2722.                         +(tda(i,j,k-1,ih,im)*rr(i,j,k,im)   &
    2723. 

  2723.                         *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
    2724. 

  2724.           rsa(i,j,k,ih,im)= rs(i,j,k,im)+ts(i,j,k,im)  &
    2725. 

  2725.                         *rsa(i,j,k-1,ih,im)*denm
    2726. 

  2726.          enddo
    2727. 

  2727.         enddo
    2728. 

  2728.        enddo
    2729. 

  2729. 

  2730.   10  continue
    2731. 

  2731. 

  2732. !-----layers are added one at a time, going up from the surface.
    2733. 

  2733. !     rra is the composite reflectance illuminated by beam radiation
    2734. 

  2734. !     rxa is the composite reflectance illuminated from above
    2735. 

  2735. !         by diffuse radiation
    2736. 

  2736. !     rra and rxa are computed from eqs. (4a) and (3a) of Chou
    2737. 

  2737.  
    2738. 

  2738. !-----for the low clouds
    2739. 

  2739. 

  2740.       do 20 is=itm,2
    2741. 

  2741. 

  2742.        do j= 1, n
    2743. 

  2743.         do i= 1, m
    2744. 

  2744.          rra(i,j,np+1,1,is)=rr(i,j,np+1,is)
    2745. 

  2745.          rxa(i,j,np+1,1,is)=rs(i,j,np+1,is)
    2746. 

  2746.          rra(i,j,np+1,2,is)=rr(i,j,np+1,is)
    2747. 

  2747.          rxa(i,j,np+1,2,is)=rs(i,j,np+1,is)
    2748. 

  2748.         enddo
    2749. 

  2749.        enddo
    2750. 

  2750. 

  2751.        do j= 1, n
    2752. 

  2752.         do i= 1, m
    2753. 

  2753.          do k=np,icb(i,j),-1
    2754. 

  2754.           denm=ts(i,j,k,is)/( 1.-rs(i,j,k,is)*rxa(i,j,k+1,1,is) )
    2755. 

  2755.           rra(i,j,k,1,is)=rr(i,j,k,is)+(td(i,j,k,is) &
    2756. 

  2756.               *rra(i,j,k+1,1,is)+tt(i,j,k,is)*rxa(i,j,k+1,1,is))*denm
    2757. 

  2757.           rxa(i,j,k,1,is)= rs(i,j,k,is)+ts(i,j,k,is) &
    2758. 

  2758.               *rxa(i,j,k+1,1,is)*denm
    2759. 

  2759.           rra(i,j,k,2,is)=rra(i,j,k,1,is)
    2760. 

  2760.           rxa(i,j,k,2,is)=rxa(i,j,k,1,is)
    2761. 

  2761.          enddo
    2762. 

  2762.         enddo
    2763. 

  2763.        enddo
    2764. 

  2764. 

  2765. !-----for middle clouds
    2766. 

  2766. 

  2767.       do 20 im=itm,2
    2768. 

  2768. 

  2769.        do j= 1, n
    2770. 

  2770.         do i= 1, m
    2771. 

  2771.          do k= icb(i,j)-1,ict(i,j),-1
    2772. 

  2772.           denm=ts(i,j,k,im)/( 1.-rs(i,j,k,im)*rxa(i,j,k+1,im,is) )
    2773. 

  2773.           rra(i,j,k,im,is)= rr(i,j,k,im)+(td(i,j,k,im)  &
    2774. 

  2774.               *rra(i,j,k+1,im,is)+tt(i,j,k,im)*rxa(i,j,k+1,im,is))*denm
    2775. 

  2775.           rxa(i,j,k,im,is)= rs(i,j,k,im)+ts(i,j,k,im)   &
    2776. 

  2776.               *rxa(i,j,k+1,im,is)*denm
    2777. 

  2777.          enddo
    2778. 

  2778.         enddo
    2779. 

  2779.        enddo
    2780. 

  2780. 

  2781.   20  continue
    2782. 

  2782. 

  2783. !-----integration over eight sky situations.
    2784. 

  2784. !     ih, im, is denotes high, middle and low cloud groups.
    2785. 

  2785. 

  2786.       do 100 ih=itm,2
    2787. 

  2787. 

  2788. !-----clear portion 
    2789. 

  2789. 

  2790.          if(ih.eq.1) then
    2791. 

  2791.            do j=1,n
    2792. 

  2792.             do i=1,m
    2793. 

  2793.              ch(i,j)=1.0-cc(i,j,1)
    2794. 

  2794.             enddo
    2795. 

  2795.            enddo
    2796. 

  2796. 

  2797.           else
    2798. 

  2798. 

  2799. !-----cloudy portion
    2800. 

  2800. 

  2801.            do j=1,n
    2802. 

  2802.             do i=1,m
    2803. 

  2803.              ch(i,j)=cc(i,j,1)
    2804. 

  2804.             enddo
    2805. 

  2805.            enddo
    2806. 

  2806. 

  2807.           endif
    2808. 

  2808. 

  2809.       do 100 im=itm,2
    2810. 

  2810. 

  2811. !-----clear portion
    2812. 

  2812. 

  2813.          if(im.eq.1) then
    2814. 

  2814. 

  2815.            do j=1,n
    2816. 

  2816.             do i=1,m
    2817. 

  2817.               cm(i,j)=ch(i,j)*(1.0-cc(i,j,2))
    2818. 

  2818.             enddo
    2819. 

  2819.            enddo
    2820. 

  2820. 

  2821.          else
    2822. 

  2822. 

  2823. !-----cloudy portion
    2824. 

  2824. 

  2825.            do j=1,n
    2826. 

  2826.             do i=1,m
    2827. 

  2827.               cm(i,j)=ch(i,j)*cc(i,j,2) 
    2828. 

  2828.             enddo
    2829. 

  2829.            enddo
    2830. 

  2830. 

  2831.          endif
    2832. 

  2832. 

  2833.       do 100 is=itm,2
    2834. 

  2834. 

  2835. !-----clear portion
    2836. 

  2836. 

  2837.          if(is.eq.1) then
    2838. 

  2838. 

  2839.            do j=1,n
    2840. 

  2840.             do i=1,m
    2841. 

  2841.              ct(i,j)=cm(i,j)*(1.0-cc(i,j,3)) 
    2842. 

  2842.             enddo
    2843. 

  2843.            enddo
    2844. 

  2844. 

  2845.          else
    2846. 

  2846. 

  2847. !-----cloudy portion
    2848. 

  2848. 

  2849.            do j=1,n
    2850. 

  2850.             do i=1,m
    2851. 

  2851.              ct(i,j)=cm(i,j)*cc(i,j,3)
    2852. 

  2852.             enddo
    2853. 

  2853.            enddo
    2854. 

  2854. 

  2855.          endif
    2856. 

  2856. 

  2857. !-----add one layer at a time, going down.
    2858. 

  2858. 

  2859.        do j= 1, n
    2860. 

  2860.         do i= 1, m
    2861. 

  2861.          do k= icb(i,j), np
    2862. 

  2862.           denm = ts(i,j,k,is)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,is) )
    2863. 

  2863.           tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,is)
    2864. 

  2864.           tta(i,j,k,ih,im)=  tda(i,j,k-1,ih,im)*tt(i,j,k,is) &
    2865. 

  2865.                +(tda(i,j,k-1,ih,im)*rr(i,j,k,is) &
    2866. 

  2866.                *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
    2867. 

  2867.           rsa(i,j,k,ih,im)= rs(i,j,k,is)+ts(i,j,k,is) &
    2868. 

  2868.                *rsa(i,j,k-1,ih,im)*denm
    2869. 

  2869.          enddo
    2870. 

  2870.         enddo
    2871. 

  2871.        enddo
    2872. 

  2872. 

  2873. !-----add one layer at a time, going up.
    2874. 

  2874. 

  2875.        do j= 1, n
    2876. 

  2876.         do i= 1, m
    2877. 

  2877.          do k= ict(i,j)-1,1,-1
    2878. 

  2878.           denm =ts(i,j,k,ih)/(1.-rs(i,j,k,ih)*rxa(i,j,k+1,im,is))
    2879. 

  2879.           rra(i,j,k,im,is)= rr(i,j,k,ih)+(td(i,j,k,ih)  &
    2880. 

  2880.               *rra(i,j,k+1,im,is)+tt(i,j,k,ih)*rxa(i,j,k+1,im,is))*denm
    2881. 

  2881.           rxa(i,j,k,im,is)= rs(i,j,k,ih)+ts(i,j,k,ih)   &
    2882. 

  2882.               *rxa(i,j,k+1,im,is)*denm
    2883. 

  2883.          enddo
    2884. 

  2884.         enddo
    2885. 

  2885.        enddo
    2886. 

  2886. 

  2887. !-----compute fluxes following eq (5) of Chou (1992)
    2888. 

  2888.  
    2889. 

  2889. !     fdndir is the direct  downward flux
    2890. 

  2890. !     fdndif is the diffuse downward flux
    2891. 

  2891. !     fupdif is the diffuse upward flux
    2892. 

  2892. 

  2893.       do k=2,np+1
    2894. 

  2894.        do j=1, n
    2895. 

  2895.         do i=1, m
    2896. 

  2896.          denm= 1./(1.- rxa(i,j,k,im,is)*rsa(i,j,k-1,ih,im))
    2897. 

  2897.          fdndir(i,j)= tda(i,j,k-1,ih,im)
    2898. 

  2898.          xx = tda(i,j,k-1,ih,im)*rra(i,j,k,im,is)
    2899. 

  2899.          fdndif(i,j)= (xx*rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm
    2900. 

  2900.          fupdif= (xx+tta(i,j,k-1,ih,im)*rxa(i,j,k,im,is))*denm
    2901. 

  2901.          flxdn(i,j,k)=fdndir(i,j)+fdndif(i,j)-fupdif
    2902. 

  2902.          flxdnu(i,j,k)=-fupdif
    2903. 

  2903.          flxdnd(i,j,k)=fdndir(i,j)+fdndif(i,j)
    2904. 

  2904.         enddo
    2905. 

  2905.        enddo
    2906. 

  2906.       enddo
    2907. 

  2907. 

  2908.        do j=1, n
    2909. 

  2909.         do i=1, m
    2910. 

  2910.          flxdn(i,j,1)=1.0-rra(i,j,1,im,is)
    2911. 

  2911.          flxdnu(i,j,1)=-rra(i,j,1,im,is)
    2912. 

  2912.          flxdnd(i,j,1)=1.0
    2913. 

  2913.         enddo
    2914. 

  2914.        enddo
    2915. 

  2915. 

  2916. !-----summation of fluxes over all (eight) sky situations.
    2917. 

  2917. 

  2918.        do k=1,np+1
    2919. 

  2919.         do j=1,n
    2920. 

  2920.          do i=1,m
    2921. 

  2921.            if(ih.eq.1 .and. im.eq.1 .and. is.eq.1) then
    2922. 

  2922.              fclr(i,j,k)=flxdn(i,j,k)
    2923. 

  2923.            endif
    2924. 

  2924.              fall(i,j,k)=fall(i,j,k)+flxdn(i,j,k)*ct(i,j)
    2925. 

  2925.              fallu(i,j,k)=fallu(i,j,k)+flxdnu(i,j,k)*ct(i,j)
    2926. 

  2926.              falld(i,j,k)=falld(i,j,k)+flxdnd(i,j,k)*ct(i,j)
    2927. 

  2927.          enddo
    2928. 

  2928.         enddo
    2929. 

  2929.        enddo
    2930. 

  2930. 

  2931.         do j=1,n
    2932. 

  2932.          do i=1,m
    2933. 

  2933.             fsdir(i,j)=fsdir(i,j)+fdndir(i,j)*ct(i,j)
    2934. 

  2934.             fsdif(i,j)=fsdif(i,j)+fdndif(i,j)*ct(i,j)
    2935. 

  2935.          enddo
    2936. 

  2936.         enddo
    2937. 

  2937. 

  2938.   100 continue
    2939. 

  2939. 

  2940.       end subroutine cldflx
    2941. 

  2941. 

  2942. !*****************************************************************
    2943. 

  2943. 

  2944.       subroutine flxco2(m,n,np,swc,swh,csm,df)
    2945. 

  2945. 

  2946. !*****************************************************************
    2947. 

  2947. 

  2948. !-----compute the reduction of clear-sky downward solar flux
    2949. 

  2949. !     due to co2 absorption.
    2950. 

  2950. 

  2951.       implicit none
    2952. 

  2952. 

  2953. !-----input parameters
    2954. 

  2954. 

  2955.       integer m,n,np
    2956. 

  2956.       real csm(m,n),swc(m,n,np+1),swh(m,n,np+1),cah(22,19)
    2957. 

  2957. 

  2958. !-----output (undated) parameter
    2959. 

  2959. 

  2960.       real df(m,n,np+1)
    2961. 

  2961. 

  2962. !-----temporary array
    2963. 

  2963. 

  2964.       integer i,j,k,ic,iw 
    2965. 

  2965.       real xx,clog,wlog,dc,dw,x1,x2,y2
    2966. 

  2966. 

  2967. !********************************************************************
    2968. 

  2968. !-----include co2 look-up table
    2969. 

  2969. 

  2970.       data ((cah(i,j),i=1,22),j= 1, 5)/ &                                     
    2971. 

  2971.        0.9923, 0.9922, 0.9921, 0.9920, 0.9916, 0.9910, 0.9899, 0.9882, &      
    2972. 

  2972.        0.9856, 0.9818, 0.9761, 0.9678, 0.9558, 0.9395, 0.9188, 0.8945, &      
    2973. 

  2973.        0.8675, 0.8376, 0.8029, 0.7621, 0.7154, 0.6647, 0.9876, 0.9876, &      
    2974. 

  2974.        0.9875, 0.9873, 0.9870, 0.9864, 0.9854, 0.9837, 0.9811, 0.9773, &      
    2975. 

  2975.        0.9718, 0.9636, 0.9518, 0.9358, 0.9153, 0.8913, 0.8647, 0.8350, &      
    2976. 

  2976.        0.8005, 0.7599, 0.7133, 0.6627, 0.9808, 0.9807, 0.9806, 0.9805, &      
    2977. 

  2977.        0.9802, 0.9796, 0.9786, 0.9769, 0.9744, 0.9707, 0.9653, 0.9573, &      
    2978. 

  2978.        0.9459, 0.9302, 0.9102, 0.8866, 0.8604, 0.8311, 0.7969, 0.7565, &      
    2979. 

  2979.        0.7101, 0.6596, 0.9708, 0.9708, 0.9707, 0.9705, 0.9702, 0.9697, &      
    2980. 

  2980.        0.9687, 0.9671, 0.9647, 0.9612, 0.9560, 0.9483, 0.9372, 0.9221, &      
    2981. 

  2981.        0.9027, 0.8798, 0.8542, 0.8253, 0.7916, 0.7515, 0.7054, 0.6551, &      
    2982. 

  2982.        0.9568, 0.9568, 0.9567, 0.9565, 0.9562, 0.9557, 0.9548, 0.9533, &      
    2983. 

  2983.        0.9510, 0.9477, 0.9428, 0.9355, 0.9250, 0.9106, 0.8921, 0.8700, &      
    2984. 

  2984.        0.8452, 0.8171, 0.7839, 0.7443, 0.6986, 0.6486/                        
    2985. 

  2985.  
    2986. 

  2986.       data ((cah(i,j),i=1,22),j= 6,10)/  &                                    
    2987. 

  2987.        0.9377, 0.9377, 0.9376, 0.9375, 0.9372, 0.9367, 0.9359, 0.9345, &      
    2988. 

  2988.        0.9324, 0.9294, 0.9248, 0.9181, 0.9083, 0.8948, 0.8774, 0.8565, &      
    2989. 

  2989.        0.8328, 0.8055, 0.7731, 0.7342, 0.6890, 0.6395, 0.9126, 0.9126, &      
    2990. 

  2990.        0.9125, 0.9124, 0.9121, 0.9117, 0.9110, 0.9098, 0.9079, 0.9052, &      
    2991. 

  2991.        0.9012, 0.8951, 0.8862, 0.8739, 0.8579, 0.8385, 0.8161, 0.7900, &      
    2992. 

  2992.        0.7585, 0.7205, 0.6760, 0.6270, 0.8809, 0.8809, 0.8808, 0.8807, &      
    2993. 

  2993.        0.8805, 0.8802, 0.8796, 0.8786, 0.8770, 0.8747, 0.8712, 0.8659, &      
    2994. 

  2994.        0.8582, 0.8473, 0.8329, 0.8153, 0.7945, 0.7697, 0.7394, 0.7024, &      
    2995. 

  2995.        0.6588, 0.6105, 0.8427, 0.8427, 0.8427, 0.8426, 0.8424, 0.8422, &      
    2996. 

  2996.        0.8417, 0.8409, 0.8397, 0.8378, 0.8350, 0.8306, 0.8241, 0.8148, &      
    2997. 

  2997.        0.8023, 0.7866, 0.7676, 0.7444, 0.7154, 0.6796, 0.6370, 0.5897, &      
    2998. 

  2998.        0.7990, 0.7990, 0.7990, 0.7989, 0.7988, 0.7987, 0.7983, 0.7978, &      
    2999. 

  2999.        0.7969, 0.7955, 0.7933, 0.7899, 0.7846, 0.7769, 0.7664, 0.7528, &      
    3000. 

  3000.        0.7357, 0.7141, 0.6866, 0.6520, 0.6108, 0.5646/                        
    3001. 

  3001.  
    3002. 

  3002.       data ((cah(i,j),i=1,22),j=11,15)/  &                                    
    3003. 

  3003.        0.7515, 0.7515, 0.7515, 0.7515, 0.7514, 0.7513, 0.7511, 0.7507, &      
    3004. 

  3004.        0.7501, 0.7491, 0.7476, 0.7450, 0.7409, 0.7347, 0.7261, 0.7144, &      
    3005. 

  3005.        0.6992, 0.6793, 0.6533, 0.6203, 0.5805, 0.5357, 0.7020, 0.7020, &      
    3006. 

  3006.        0.7020, 0.7019, 0.7019, 0.7018, 0.7017, 0.7015, 0.7011, 0.7005, &      
    3007. 

  3007.        0.6993, 0.6974, 0.6943, 0.6894, 0.6823, 0.6723, 0.6588, 0.6406, &      
    3008. 

  3008.        0.6161, 0.5847, 0.5466, 0.5034, 0.6518, 0.6518, 0.6518, 0.6518, &      
    3009. 

  3009.        0.6518, 0.6517, 0.6517, 0.6515, 0.6513, 0.6508, 0.6500, 0.6485, &      
    3010. 

  3010.        0.6459, 0.6419, 0.6359, 0.6273, 0.6151, 0.5983, 0.5755, 0.5458, &      
    3011. 

  3011.        0.5095, 0.4681, 0.6017, 0.6017, 0.6017, 0.6017, 0.6016, 0.6016, &      
    3012. 

  3012.        0.6016, 0.6015, 0.6013, 0.6009, 0.6002, 0.5989, 0.5967, 0.5932, &      
    3013. 

  3013.        0.5879, 0.5801, 0.5691, 0.5535, 0.5322, 0.5043, 0.4700, 0.4308, &      
    3014. 

  3014.        0.5518, 0.5518, 0.5518, 0.5518, 0.5518, 0.5518, 0.5517, 0.5516, &      
    3015. 

  3015.        0.5514, 0.5511, 0.5505, 0.5493, 0.5473, 0.5441, 0.5393, 0.5322, &      
    3016. 

  3016.        0.5220, 0.5076, 0.4878, 0.4617, 0.4297, 0.3929/                        
    3017. 

  3017.  
    3018. 

  3018.       data ((cah(i,j),i=1,22),j=16,19)/ &                                     
    3019. 

  3019.        0.5031, 0.5031, 0.5031, 0.5031, 0.5031, 0.5030, 0.5030, 0.5029, &      
    3020. 

  3020.        0.5028, 0.5025, 0.5019, 0.5008, 0.4990, 0.4960, 0.4916, 0.4850, &      
    3021. 

  3021.        0.4757, 0.4624, 0.4441, 0.4201, 0.3904, 0.3564, 0.4565, 0.4565, &      
    3022. 

  3022.        0.4565, 0.4564, 0.4564, 0.4564, 0.4564, 0.4563, 0.4562, 0.4559, &      
    3023. 

  3023.        0.4553, 0.4544, 0.4527, 0.4500, 0.4460, 0.4400, 0.4315, 0.4194, &      
    3024. 

  3024.        0.4028, 0.3809, 0.3538, 0.3227, 0.4122, 0.4122, 0.4122, 0.4122, &      
    3025. 

  3025.        0.4122, 0.4122, 0.4122, 0.4121, 0.4120, 0.4117, 0.4112, 0.4104, &      
    3026. 

  3026.        0.4089, 0.4065, 0.4029, 0.3976, 0.3900, 0.3792, 0.3643, 0.3447, &      
    3027. 

  3027.        0.3203, 0.2923, 0.3696, 0.3696, 0.3696, 0.3696, 0.3696, 0.3696, &      
    3028. 

  3028.        0.3695, 0.3695, 0.3694, 0.3691, 0.3687, 0.3680, 0.3667, 0.3647, &      
    3029. 

  3029.        0.3615, 0.3570, 0.3504, 0.3409, 0.3279, 0.3106, 0.2892, 0.2642/        
    3030. 

  3030. 

  3031. !********************************************************************
    3032. 

  3032. !-----table look-up for the reduction of clear-sky solar
    3033. 

  3033. !     radiation due to co2. The fraction 0.0343 is the
    3034. 

  3034. !     extraterrestrial solar flux in the co2 bands.
    3035. 

  3035. 

  3036.       do k= 2, np+1
    3037. 

  3037.        do j= 1, n
    3038. 

  3038.         do i= 1, m
    3039. 

  3039.           xx=1./.3
    3040. 

  3040.           clog=log10(swc(i,j,k)*csm(i,j))
    3041. 

  3041.           wlog=log10(swh(i,j,k)*csm(i,j))
    3042. 

  3042.           ic=int( (clog+3.15)*xx+1.)
    3043. 

  3043.           iw=int( (wlog+4.15)*xx+1.)
    3044. 

  3044.           if(ic.lt.2)ic=2
    3045. 

  3045.           if(iw.lt.2)iw=2
    3046. 

  3046.           if(ic.gt.22)ic=22
    3047. 

  3047.           if(iw.gt.19)iw=19     
    3048. 

  3048.           dc=clog-float(ic-2)*.3+3.
    3049. 

  3049.           dw=wlog-float(iw-2)*.3+4.   
    3050. 

  3050.           x1=cah(1,iw-1)+(cah(1,iw)-cah(1,iw-1))*xx*dw
    3051. 

  3051.           x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw
    3052. 

  3052.           y2=x2+(cah(ic,iw-1)-cah(ic-1,iw-1))*xx*dc
    3053. 

  3053.           if (x1.lt.y2) x1=y2
    3054. 

  3054.           df(i,j,k)=df(i,j,k)+0.0343*(x1-y2)
    3055. 

  3055.         enddo     
    3056. 

  3056.        enddo      
    3057. 

  3057.       enddo      
    3058. 

  3058. 

  3059.       end subroutine flxco2
    3060. 

  3060. 

  3061. !*****************************************************************
    3062. 

  3062. 

  3063.       subroutine o3prof (np, pres, ozone, its, ite, kts, kte, p, o3)
    3064. 

  3064. 

  3065. !*****************************************************************
    3066. 

  3066.       implicit none
    3067. 

  3067. !*****************************************************************
    3068. 

  3068. !
    3069. 

  3069.       integer iprof,m,np,its,ite,kts,kte
    3070. 

  3070.       integer i,k,ko,kk
    3071. 

  3071.       real pres(np),ozone(np)
    3072. 

  3072.       real p(its:ite,kts:kte),o3(its:ite,kts:kte)
    3073. 

  3073.  
    3074. 

  3074. !     Statement function 
    3075. 

  3075.  
    3076. 

  3076.       real Linear, x1, y1, x2, y2, x
    3077. 

  3077.       Linear(x1, y1, x2, y2, x) =  &
    3078. 

  3078.             (y1 * (x2 - x) + y2 * (x - x1)) / (x2 - x1)
    3079. 

  3079. !
    3080. 

  3080.       do k = 1,np
    3081. 

  3081.         pres(k) = alog(pres(k))
    3082. 

  3082.       enddo
    3083. 

  3083.       do k = kts,kte
    3084. 

  3084.         do i = its, ite
    3085. 

  3085.           p(i,k) = alog(p(i,k))
    3086. 

  3086.         end do
    3087. 

  3087.       end do
    3088. 

  3088. 

  3089. ! assume the pressure at model top is greater than pres(1)
    3090. 

  3090. ! if it is not, this part needs to change
    3091. 

  3091. 

  3092.       do i = its, ite
    3093. 

  3093.         ko = 1
    3094. 

  3094.         do k = kts+1, kte
    3095. 

  3095.           do while (ko .lt. np .and. p(i,k) .gt. pres(ko))
    3096. 

  3096.             ko = ko + 1
    3097. 

  3097.           end do
    3098. 

  3098.           o3(i,k) =  Linear (pres(ko),   ozone(ko),    &
    3099. 

  3099.                              pres(ko-1), ozone(ko-1),  &
    3100. 

  3100.                              p(i,k))
    3101. 

  3101.           ko = ko - 1
    3102. 

  3102.         end do
    3103. 

  3103.       end do
    3104. 

  3104. 

  3105. ! calculate top lay O3
    3106. 

  3106. 

  3107.       do i = its, ite
    3108. 

  3108.         ko = 1
    3109. 

  3109.         k = kts
    3110. 

  3110.         do while (ko .le. np .and. p(i,k) .gt. pres(ko))
    3111. 

  3111.            ko = ko + 1
    3112. 

  3112.         end do
    3113. 

  3113.         IF (ko-1 .le. 1) then
    3114. 

  3114.            O3(i,k)=ozone(k)
    3115. 

  3115.         ELSE
    3116. 

  3116.            O3(i,k)=0.
    3117. 

  3117.            do kk=ko-2,1,-1
    3118. 

  3118.               O3(i,k)=O3(i,k)+ozone(kk)*(pres(kk+1)-pres(kk))
    3119. 

  3119.            enddo
    3120. 

  3120.            O3(i,k)=O3(i,k)/(pres(ko-1)-pres(1))
    3121. 

  3121.         ENDIF
    3122. 

  3122. !       print*,'O3=',i,k,ko,O3(i,k),p(i,k),ko,pres(ko),pres(ko-1)
    3123. 

  3123.       end do
    3124. 

  3124.       
    3125. 

  3125.       end subroutine o3prof
    3126. 

  3126. 

  3127. !-----------------------------------------
    3128. 

  3128.     SUBROUTINE gsfc_swinit(cen_lat, allowed_to_read)
    3129. 

  3129. 

  3130.     REAL, INTENT(IN    )      ::        cen_lat
    3131. 

  3131.     LOGICAL, INTENT(IN    )   ::       allowed_to_read
    3132. 

  3132. 

  3133.         center_lat=cen_lat
    3134. 

  3134. 

  3135.     END SUBROUTINE gsfc_swinit
    3136. 

  3136. 

  3137. 

  3138. END MODULE module_ra_gsfcsw
    3139. 

  3139.