File: C:\NOAA\NEMS_11731\src\atmos\phys\module_RA_RRTM.F90

1     
2     !
3           MODULE MODULE_RA_RRTM
4     !
5     !-----------------------------------------------------------------------
6     !
7     !***  THE RADIATION DRIVERS AND PACKAGES
8     !
9     !-----------------------------------------------------------------------
10     !
11           USE MODULE_INCLUDE
12     !
13           USE MODULE_CONSTANTS,ONLY : R,CP,PI,EPSQ,STBOLT,EP_2
14           USE MODULE_n_RADIATION_DRIVER,   ONLY : n_RADINIT, n_GRRAD
15           USE MODULE_n_RADIATION_ASTRONOMY,ONLY : n_ASTRONOMY
16     !     USE n_OZNE_DEF, ONLY: LEVOZC, LATSOZP, BLATC, TIMEOZC, TIMEOZ,      &
17           USE OZNE_DEF, ONLY: LEVOZC, LATSOZP, BLATC, TIMEOZC, TIMEOZ,      &
18                               KOZC, DPHIOZC, LATSOZC, PL_COEFF, LEVOZP
19     
20           USE MODULE_MP_ETANEW,ONLY : RHgrd,T_ICE,FPVS
21     
22     !-----------------------------------------------------------------------
23     !
24           IMPLICIT NONE
25     !
26     !-----------------------------------------------------------------------
27     !
28           PRIVATE
29     !
30           PUBLIC :: RRTM,RRTM_INIT
31     !
32     !-----------------------------------------------------------------------
33     !
34     !--- Used for Gaussian look up tables
35     !
36           REAL, PRIVATE,PARAMETER :: XSDmax=3.1, DXSD=.01
37           INTEGER, PRIVATE,PARAMETER :: NXSD=XSDmax/DXSD
38           REAL, DIMENSION(NXSD),PRIVATE,SAVE :: AXSD
39           REAL, PRIVATE :: RSQR
40           LOGICAL, PRIVATE,SAVE :: SDprint=.FALSE.
41     
42     !-------------------------------
43           INTEGER, SAVE, DIMENSION(3)     :: LTOP
44           REAL,SAVE,DIMENSION(4) :: PTOPC
45     !--------------------------------
46     !
47           REAL, PARAMETER ::  &
48          &   TRAD_ice=0.5*T_ice      & !--- Very tunable parameter
49          &,  ABSCOEF_W=800.            & !--- Very tunable parameter
50          &,  ABSCOEF_I=0.            & !--- Very tunable parameter
51          &,  Qconv=0.01e-3            & !--- Very tunable parameter
52     
53          &,  CTauCW=ABSCOEF_W*Qconv  &
54          &,  CTauCI=ABSCOEF_I*Qconv
55     
56     
57     
58     !-----------------------------------------------------------------------
59     !-----------------------------------------------------------------------
60     !***  THE RADIATION PACKAGE OPTIONS
61     !-----------------------------------------------------------------------
62     !
63           CONTAINS
64     !
65     !-----------------------------------------------------------------------
66     !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
67     !-----------------------------------------------------------------------
68           SUBROUTINE RRTM (NTIMESTEP,DT_INT,JDAT                            &
69          &                    ,NPHS,GLAT,GLON                               &
70          &                    ,NRADS,NRADL                                  &
71          &                    ,DSG2,SGML2,PDSG1,PSGML1                      &
72          &                    ,PT,PD                                        &
73          &                    ,T,Q,CW,O3                                    &
74          &                    ,ALBEDO                                       &
75          &                    ,F_ICE,F_RAIN                                 &
76          &                    ,P_QV,P_QC,P_QR,P_QI,P_QS,P_QG                &
77          &                    ,SM,CLDFRA                                    &
78          &                    ,NUM_WATER,WATER                              &
79          &                    ,RLWTT,RSWTT                                  &
80          &                    ,RLWIN,RSWIN                                  &
81          &                    ,RSWINC,RSWOUT                                &
82          &                    ,RLWTOA,RSWTOA                                &
83          &                    ,CZMEAN,SIGT4                                 &
84          &                    ,CFRACL,CFRACM,CFRACH                         &
85          &                    ,ACFRST,NCFRST                                &
86          &                    ,ACFRCV,NCFRCV                                &
87          &                    ,CUPPT,SNOW,SI                                &
88          &                    ,HTOP,HBOT                                    &
89          &                    ,TSKIN,Z0,SICE,F_RIMEF,MXSNAL,SGM,STDH,OMGALF &
90          &                    ,IMS,IME,JMS,JME                              &
91          &                    ,ITS,ITE,JTS,JTE                              &
92          &                    ,LM                                           &
93          &                    ,MYPE )
94     !
95     !-----------------------------------------------------------------------
96     !
97           IMPLICIT NONE
98     !
99     !-----------------------------------------------------------------------
100     !
101           INTEGER,INTENT(IN) :: IME,IMS,ITE,ITS                             &
102          &                     ,JME,JMS,JTE,JTS                             &
103          &                     ,LM,MYPE                                     &
104          &                     ,NTIMESTEP,DT_INT                            &
105          &                     ,NPHS,NRADL,NRADS                            &
106          &                     ,NUM_WATER                      
107     !
108           INTEGER,INTENT(IN) :: JDAT(8)
109     !
110           INTEGER,INTENT(IN) :: P_QV,P_QC,P_QR,P_QI,P_QS,P_QG
111     !
112           INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: NCFRCV,NCFRST
113     !
114           REAL,INTENT(IN) :: PT
115     !
116           REAL,DIMENSION(1:LM),INTENT(IN) :: DSG2,PDSG1,PSGML1,SGML2
117     !
118           REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CUPPT               &
119                                                        ,GLAT,GLON           &
120                                                        ,PD,SM,SNOW,SI        
121     !
122           REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: ALBEDO
123     
124           REAL,DIMENSION(IMS:IME,JMS:JME,LM),INTENT(IN) :: CW,O3,Q,T
125     !
126           REAL,DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(IN) :: F_ICE,F_RAIN
127     !
128           REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: ACFRCV,ACFRST    &
129                                                           ,RLWIN,RLWTOA     &
130                                                           ,RSWIN,RSWOUT     &
131                                                           ,HBOT,HTOP        &
132                                                           ,RSWINC,RSWTOA
133     !
134           REAL,DIMENSION(IMS:IME,JMS:JME,LM),INTENT(INOUT) :: RLWTT,RSWTT
135     !
136           REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: CFRACH,CFRACL    &
137                                                           ,CFRACM,CZMEAN    &
138                                                           ,SIGT4
139     !
140           REAL,DIMENSION(IMS:IME,JMS:JME,LM,NUM_WATER),INTENT(INOUT) :: WATER
141     !
142           REAL,DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(OUT) :: CLDFRA
143     !
144            REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: TSKIN,Z0,SICE      &
145                                                         ,MXSNAL,STDH        
146     !
147            REAL,DIMENSION(1:LM+1),INTENT(IN) :: SGM
148     !
149            REAL,DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(IN) :: F_RIMEF,OMGALF
150     !
151     !-----------------------------------------------------------------------
152     !***  LOCAL VARIABLES
153     !-----------------------------------------------------------------------
154     !
155           LOGICAL :: LSLWR, LSSAV, LDIAG3D, LPRNT, SASHAL, LSSWR
156     !
157           INTEGER,PARAMETER :: IFLIP=0                                         
158     !
159           INTEGER :: NP3D, ISOL, ICO2, ICWP, IALB, IEMS, IAER, ALBTYPE, LONR, &
160                      LATS_NODE_R, LATR, IPT_LATS_NODE_R, NFXR, NTRAC, KFLIP,  &
161                      I, L, J, K, NTOZ, NCLDX, NTCW, IOVR_SW, IOVR_LW
162     !
163           INTEGER,SAVE :: K1OZ, K2OZ
164     !
165           INTEGER,DIMENSION(JTS:JTE) :: LONSPERLAR, GLOBAL_LATS_R
166     !
167           REAL*8 :: FHSWR, FHLWR, FHAER, SOLCON,SLAG, SDEC, CDEC,DTSW, DTLW, RTvR
168     !
169           REAL*8,SAVE :: FACOZ
170     !
171           REAL*8,DIMENSION(1) :: FLGMIN_L, CV, CVB, CVT,XLAT, HPRIME_V, TSEA, &
172                                  TISFC, FICE, ZORL, SLMSK, SNWDPH, SNCOVR,    &
173                                  SNOALB, ALVSF1, ALNSF1, ALVWF1, ALNWF1,      &
174                                  FACSF1, FACWF1, SFCNSW, SFCDSW, SFALB,       &
175                                  SFCDLW, TSFLW, TOAUSW, TOADSW, SFCCDSW,      &
176                                  TOAULW, SFCUSW
177     !
178           REAL*8,DIMENSION(LM) :: CLDCOV_V,PRSL,PRSLK,GT,GQ, VVEL,F_ICEC,     &
179                                    F_RAINC,R_RIME,TAUCLOUDS,CLDF
180     
181           REAL*8,DIMENSION(LM+1) :: PRSI , RSGM
182     !
183           REAL*8,DIMENSION(JTS:JTE) :: COSLAT_R, SINLAT_R
184     !
185           REAL*8,DIMENSION(ITS:ITE,JTS:JTE) :: XLON, COSZEN,      &
186                                                            COSZDG, SFCALBEDO
187     !
188           REAL*8,DIMENSION(27) :: FLUXR_V
189     !
190           REAL*8,DIMENSION(1,LM,3) :: GR1   
191     !
192           REAL*8,DIMENSION(LM) :: SWH, HLW
193     !
194           REAL,DIMENSION(ITS:ITE,JTS:JTE,27) :: FLUXR_VIJ
195           
196           REAL :: BLATC4
197     !
198           REAL,DIMENSION(ITS:ITE,JTS:JTE,1:LM+1) :: P8W
199     !
200           REAL,DIMENSION(ITS:ITE,JTS:JTE,1:LM)   :: P_PHY
201     !
202           LOGICAL,SAVE :: FIRST
203           DATA FIRST / .TRUE. /
204     !
205           INTEGER :: DAYS(13), IDAY, IMON, MIDMON, ID
206     !
207           INTEGER, SAVE :: MIDM, MIDP
208     !
209           LOGICAL :: CHANGE
210     !
211           DATA DAYS / 31,28,31,30,31,30,31,31,30,31,30,31,30 /
212     !
213           REAL :: ZEN,DZEN,ALB1,ALB2
214     !
215           INTEGER :: IR,IQ,JX
216     !
217           REAL,PARAMETER :: TWENTY=20.0, &
218                             HP537=0.537, &
219                             ONE=1., &
220                             DEGRAD1=180.0/PI, &
221                             H74E1=74.0, &
222                             HAF=0.5, &
223                             HNINETY=90., &
224                             FIFTY=50., &
225                             QUARTR=0.25, &
226                             HNINE=9.0, &
227                             HP1=0.1, &
228                             H15E1=15.0
229           REAL,DIMENSION(20) :: ZA
230           REAL,DIMENSION(19) :: DZA
231           REAL,DIMENSION(21,20) :: ALBD
232           REAL,DIMENSION(21) :: TRN
233     !
234           DATA TRN/.00,.05,.10,.15,.20,.25,.30,.35,.40,.45,.50,.55,.60,.65, &
235                    .70,.75,.80,.85,.90,.95,1.00/
236           
237           DATA  ALBD/.061,.062,.072,.087,.115,.163,.235,.318,.395,.472,.542, &
238            .604,.655,.693,.719,.732,.730,.681,.581,.453,.425,.061,.062,.070, &
239            .083,.108,.145,.198,.263,.336,.415,.487,.547,.595,.631,.656,.670, &
240            .652,.602,.494,.398,.370,.061,.061,.068,.079,.098,.130,.174,.228, &
241            .290,.357,.424,.498,.556,.588,.603,.592,.556,.488,.393,.342,.325, &
242            .061,.061,.065,.073,.086,.110,.150,.192,.248,.306,.360,.407,.444, &
243            .469,.480,.474,.444,.386,.333,.301,.290,.061,.061,.065,.070,.082, &
244            .101,.131,.168,.208,.252,.295,.331,.358,.375,.385,.377,.356,.320, &
245            .288,.266,.255,.061,.061,.063,.068,.077,.092,.114,.143,.176,.210, &
246            .242,.272,.288,.296,.300,.291,.273,.252,.237,.266,.220,.061,.061, &
247            .062,.066,.072,.084,.103,.127,.151,.176,.198,.219,.236,.245,.250, &
248            .246,.235,.222,.211,.205,.200,                                    &
249                      .061,.061,.061,.065,.071,.079,.094,.113,.134,.154,.173, &
250            .185,.190,.193,.193,.190,.188,.185,.182,.180,.178,.061,.061,.061, &
251            .064,.067,.072,.083,.099,.117,.135,.150,.160,.164,.165,.164,.162, &
252            .160,.159,.158,.157,.157,.061,.061,.061,.062,.065,.068,.074,.084, &
253            .097,.111,.121,.127,.130,.131,.131,.130,.129,.127,.126,.125,.122, &
254            .061,.061,.061,.061,.062,.064,.070,.076,.085,.094,.101,.105,.107, &
255            .106,.103,.100,.097,.096,.095,.095,.095,.061,.061,.061,.060,.061, &
256            .062,.065,.070,.075,.081,.086,.089,.090,.088,.084,.080,.077,.075, &
257            .074,.074,.074,.061,.061,.060,.060,.060,.061,.063,.065,.068,.072, &
258            .076,.077,.076,.074,.071,.067,.064,.062,.061,.061,.061,.061,.061, &
259            .060,.060,.060,.060,.061,.062,.065,.068,.069,.069,.068,.065,.061, &
260            .058,.055,.054,.053,.052,.052,                                    &
261                      .061,.061,.060,.060,.060,.060,.060,.060,.062,.065,.065, &
262            .063,.060,.057,.054,.050,.047,.046,.045,.044,.044,.061,.061,.060, &
263            .060,.060,.059,.059,.059,.059,.059,.058,.055,.051,.047,.043,.039, &
264            .035,.033,.032,.031,.031,.061,.061,.060,.060,.060,.059,.059,.058, &
265            .057,.056,.054,.051,.047,.043,.039,.036,.033,.030,.028,.027,.026, &
266            .061,.061,.060,.060,.060,.059,.059,.058,.057,.055,.052,.049,.045, &
267            .040,.036,.032,.029,.027,.026,.025,.025,.061,.061,.060,.060,.060, &
268            .059,.059,.058,.056,.053,.050,.046,.042,.038,.034,.031,.028,.026, &
269            .025,.025,.025,.061,.061,.060,.060,.059,.058,.058,.057,.055,.053, &
270            .050,.046,.042,.038,.034,.030,.028,.029,.025,.025,.025/
271     !
272           DATA ZA/90.,88.,86.,84.,82.,80.,78.,76.,74.,70.,66.,62.,58.,54., &
273                   50.,40.,30.,20.,10.,0.0/
274     !
275           DATA DZA/8*2.0,6*4.0,5*10.0/
276     !
277           REAL :: ALBD0,ALVD1,ALND1
278     !
279           REAL,DIMENSION(1) :: ALVB,ALNB,ALVD,ALND
280     !
281           INTEGER :: IRTN, IERROR, o3clm_unit
282     !
283           REAL :: RJDAY, WEI2M, WEI1M, WEI1S, WEI2S, BLTO, BLNO
284     !
285           INTEGER :: JDOY, JDAY, JDOW, MMM, MMP, MM, IRET, MONEND, &
286                      MON1, IS2, ISX, KPD9, IS1, NN, MON2, MON, IS, &  
287                      LUGB, LEN, M1, M2, K1, K2, JMSK, IMSK       
288     !
289           INTEGER :: KPDALB(4)
290     !
291           CHARACTER*500 :: FNALBC,FNALBC2,FNMSKH
292     !
293           REAL :: ALBCLM(1,4), ALFCLM(1,2),  &
294                   DAYHF(13)
295     !
296           DATA DAYHF/ 15.5, 45.0, 74.5,105.0,135.5,166.0,          &
297                     196.5,227.5,258.0,288.5,319.0,349.5,380.5/
298     !
299           REAL,ALLOCATABLE :: ALB(:,:,:), ALF(:,:)
300     !
301           INTEGER :: MON1S, MON2S, SEA1S, SEA2S, SEA1, SEA2
302     
303           DATA MON1S/0/, MON2S/0/, SEA1S/0/, SEA2S/0/
304     !
305           SAVE  ALB, ALF, MON1S, MON2S, SEA1S, SEA2S, DAYHF, K1, K2, M1, M2
306     !
307           REAL :: ALBLMX,ALBLMN,ALBOMX,ALBOMN,ALBSMX, &
308                   ALBSMN,ALBIMX,ALBIMN,ALBJMX,ALBJMN, &
309                   EPSALB,PERCRIT
310     !
311           REAL :: WV,QICE,QCLD,CLFR,ESAT,QSAT,RHUM,RHtot,ARG,SDM,   &
312                    PMOD,CONVPRATE,CLSTP,P1,P2,CC1,CC2,CLDMAX,CL1,CL2, &
313                    CR1,DPCL,PRS1,PRS2,DELP,TCLD,CTau,CFSmax,CFCmax,  &
314                    CFRAVG,TDUM
315     !
316           INTEGER :: IXSD,NTSPH,NRADPP,NC,NMOD,LCNVT,LCNVB,NLVL,MALVL, &
317                      LLTOP,LLBOT,KBT2,KTH1,KBT1,KTH2,KTOP1,LM1,LL
318     !
319           REAL, PARAMETER :: EPSQ1=1.E-5,EPSQ=1.E-12,EPSO3=1.E-10,H0=0., &
320                              H1=1.,HALF=.5,T0C=273.15,CUPRATE=24.*1000., &
321                              HPINC=HALF*1.E1, CLFRmin=0.01, TAUCmax=4.161, &
322                              XSDmin=-XSDmax, DXSD1=-DXSD, STSDM=0.01, & 
323                              CVSDM=.04,DXSD2=HALF*DXSD,DXSD2N=-DXSD2,PCLDY=0.25
324     !
325           REAL,DIMENSION(10),SAVE :: CC,PPT
326           LOGICAL, SAVE :: CNCLD=.TRUE.
327     !
328           DATA CC/0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/
329           DATA PPT/0.,.14,.31,.70,1.6,3.4,7.7,17.,38.,85./
330     !
331           REAL,DIMENSION(0:LM)  :: CLDAMT
332     !
333           LOGICAL :: BITX,BITY,BITZ,BITW,BIT1,BIT2,NEW_CLOUD
334     !
335           REAL :: CTHK(3)
336           DATA CTHK/20000.0,20000.0,20000.0/
337     ! 
338           REAL,DIMENSION(ITS:ITE,JTS:JTE,3):: CLDCFR
339           INTEGER,DIMENSION(ITS:ITE,JTS:JTE,3):: MBOT,MTOP
340     
341           REAL,DIMENSION(ITS:ITE,JTS:JTE):: CUTOP,CUBOT
342           
343           REAL,DIMENSION(ITS:ITE,JTS:JTE,LM) :: TauCI  &
344                                                             ,CSMID,CCMID
345     !
346           INTEGER,DIMENSION(ITS:ITE,JTS:JTE,LM+1) :: KTOP, KBTM
347     
348           REAL,DIMENSION(ITS:ITE,JTS:JTE,LM+1) :: CAMT
349     !
350           INTEGER,DIMENSION(ITS:ITE,JTS:JTE) :: NCLDS, KCLD
351     !
352           REAL,DIMENSION(ITS:ITE,JTS:JTE,LM) :: TAUTOTAL
353     !
354           INTEGER :: NKTP,NBTM,NCLD,LML
355     !
356           REAL :: CLFR1,TauC,QSUM,DELPTOT
357     !--------------------------------------------------------------------------------------------------
358     !
359     !***THIS SUBROUTINE SELECTS AND PREPARES THE NECESSARY INPUTS FOR GRRAD (GFS RRTM DRIVER)
360     !
361     !   GRRAD IS CALLED COLUMN BY COLUMN
362     !
363     !INPUTS/OUPUTS OF GRRAD: 
364     !    INPUT VARIABLES:                                                   !
365     !      PRSI  (LM+1)    : MODEL LEVEL PRESSURE IN CB (KPA)               !
366     !      PRSL  (LM)      : MODEL LAYER MEAN PRESSURE IN CB (KPA)          !
367     !      PRSLK (LM)      : PRESSURE IN CB (KPA)                           !
368     !      GT    (LM)      : MODEL LAYER MEAN TEMPERATURE IN K              !
369     !      GQ    (LM)      : LAYER SPECIFIC HUMIDITY IN GM/GM               !
370     !      GR1   (LM,NTRAC): TRACER ARRAY (WATER, OZONE,....)               !
371     !      VVEL   (LM)      : LAYER MEAN VERTICAL VELOCITY IN CB/SEC         ! !not used
372     !      SLMSK (1)       : SEA/LAND MASK ARRAY (SEA:0,LAND:1,SEA-ICE:2)   !
373     !      XLON,XLAT       : GRID LONGITUDE/LATITUDE IN RADIANS             !
374     !      TSEA  (1)       : SURFACE TEMPERATURE IN K                       !
375     !      SNWDPH (1)       : SNOW DEPTH WATER EQUIVALENT IN MM              !
376     !      SNCOVR(1)       : SNOW COVER IN FRACTION                         !
377     !      SNOALB(1)       : MAXIMUM SNOW ALBEDO IN FRACTION                !
378     !      ZORL  (1)       : SURFACE ROUGHNESS IN CM                        !
379     !      HPRIM_V (1)       : TOPOGRAPHIC STANDARD DEVIATION IN M            !
380     !      ALVSF1 (1)       : MEAN VIS ALBEDO WITH STRONG COSZ DEPENDENCY    !
381     !      ALNSF1 (1)       : MEAN NIR ALBEDO WITH STRONG COSZ DEPENDENCY    !
382     !      ALVWF1 (1)       : MEAN VIS ALBEDO WITH WEAK COSZ DEPENDENCY      !
383     !      ALNWF1 (1)       : MEAN NIR ALBEDO WITH WEAK COSZ DEPENDENCY      !
384     !      FACSF1 (1)       : FRACTIONAL COVERAGE WITH STRONG COSZ DEPENDEN  !
385     !      FACWF1 (1)       : FRACTIONAL COVERAGE WITH WEAK COSZ DEPENDENCY  !
386     !      FICE  (1)       : ICE FRACTION OVER OPEN WATER GRID              !
387     !      TISFC (1)       : SURFACE TEMPERATURE OVER ICE FRACTION          !
388     !      SOLCON          : SOLAR CONSTANT (SUN-EARTH DISTANT ADJUSTED)    !
389     !      COSZEN          : MEAN COS OF ZENITH ANGLE OVER RAD CALL PERIOD  !
390     !      COSZDG          : DAYTIME MEAN COSZ OVER RAD CALL PERIOD         !
391     !      K1OZ,K2OZ,FACOZ : PARAMETERS FOR CLIMATOLOGICAL OZONE            !
392     !      CV    (1)       : FRACTION OF CONVECTIVE CLOUD                   ! !not used
393     !      CVT, CVB (1)    : CONVECTIVE CLOUD TOP/BOTTOM PRESSURE IN CB     ! !not used
394     !      IOVRSW/IOVRLW   : CONTROL FLAG FOR CLOUD OVERLAP (SW/LW RAD)     !
395     !                        =0 RANDOM OVERLAPPING CLOUDS                   !
396     !                        =1 MAX/RAN OVERLAPPING CLOUDS                  !
397     !      F_ICEC (LM)     : FRACTION OF CLOUD ICE  (IN FERRIER SCHEME)     !
398     !      F_RAINC(LM)     : FRACTION OF RAIN WATER (IN FERRIER SCHEME)     !
399     !      RRIME  (LM)     : MASS RATIO OF TOTAL TO UNRIMED ICE ( >= 1 )    !
400     !      FLGMIN_L(1)     : MINIMIM LARGE ICE FRACTION                     !
401     !      NP3D            : =3 BRAD FERRIER MICROPHYSICS SCHEME            ! !only stratiform clouds. optical prop. and cloud fractions are calculated in grrad
402     !                        =4 ZHAO/CARR/SUNDQVIST MICROPHYSICS SCHEME     ! !not used
403     !                        =5 NAM stratiform + convective cloud optical   ! !new - clouds optical depth and fraction calculated here and are input for grrad
404     !      NTCW            : =0 NO CLOUD CONDENSATE CALCULATED              !
405     !                        >0 ARRAY INDEX LOCATION FOR CLOUD CONDENSATE   !
406     !      NCLDX           : ONLY USED WHEN NTCW .GT. 0                     !
407     !      NTOZ            : =0 CLIMATOLOGICAL OZONE PROFILE                !
408     !                        >0 INTERACTIVE OZONE PROFILE                   ! !does not work currently
409     !      NTRAC           : DIMENSION VERIABLE FOR ARRAY GR1               !
410     !      NFXR            : SECOND DIMENSION OF INPUT/OUTPUT ARRAY FLUXR   !
411     !      DTLW, DTSW      : TIME DURATION FOR LW/SW RADIATION CALL IN SEC  !
412     !      LSSWR, LSLWR    : LOGICAL FLAGS FOR SW/LW RADIATION CALLS        !
413     !      LSSAV           : LOGICAL FLAG FOR STORE 3-D CLOUD FIELD         !
414     !      LDIAG3D         : LOGICAL FLAG FOR STORE 3-D DIAGNOSTIC FIELDS   !
415     !      LM              : VERTICAL LAYER DIMENSION                       !
416     !      IFLIP           : CONTROL FLAG FOR IN/OUT VERTICAL INDEXING      !
417     !                        =0 INDEX FROM TOA TO SURFACE                   !
418     !                        =1 INDEX FROM SURFACE TO TOA                   !
419     !      MYPE            : CONTROL FLAG FOR PARALLEL PROCESS              !
420     !      LPRNT           : CONTROL FLAG FOR DIAGNOSTIC PRINT OUT          !
421     !
422     !      LATSOZC,LEVOZC,BLATC,DPHIOZC,TIMEOZC: OZONE PARAMETERS FOR NMMB  ! !new
423     !      ALBTYPE         : OPTION FOR SELCTING NMMB ALBEDO                ! !new
424     !      TAUCLOUDS(LM)   : CLOUD OPTICAL DEPTH FROM NMMB (ferrier+bmj)    ! !new
425     !      CLDF(LM)        : CLOUD FRACTION FROM NMMB (ferrier+bmj)         ! !new
426     !                                                                       !
427     !    OUTPUT VARIABLES:                                                  !
428     !      SWH (LM)       : TOTAL SKY SW HEATING RATE IN K/SEC              !
429     !      SFCNSW(1)      : TOTAL SKY SURFACE NET SW FLUX IN W/M**2         !
430     !      SFCDSW(1)      : TOTAL SKY SURFACE DOWNWARD SW FLUX IN W/M**2    !
431     !      SFALB (1)      : MEAN SURFACE DIFFUSED ALBEDO                    !
432     !      HLW (LM)       : TOTAL SKY LW HEATING RATE IN K/SEC              !
433     !      SFCDLW(1)      : TOTAL SKY SURFACE DOWNWARD LW FLUX IN W/M**2    !
434     !      TSFLW (1)      : SURFACE AIR TEMP DURING LW CALCULATION IN K     !
435     !
436     !      TOAUSW (IM)    : TOTAL SKY TOA UPWARD SW FLUX IN W/M**2         ! !new
437     !      TOADSW (IM)    : TOTAL SKY TOA DOWNWARD SW FLUX IN W/M**2       ! !new
438     !      SFCCDSW(IM)    : CLEAR SKY SURFACE SW DOWNWARD FLUX IN W/M**2   ! !new
439     !      TOAULW (IM)    : TOTAL SKY TOA LW FLUX W/M**2                   ! !new
440     !      SFCUSW (IM)    : TOTAL SKY SURFACE SW UPWARD FLUX IN W/M**2     ! !new
441     !                                                                       !
442     !    INPUT AND OUTPUT VARIABLES:                                        !
443     !      FLUXR_V (IX,NFXR) : TO SAVE 2-D FIELDS                             !
444     !      CLDCOV_V(IX,LM)   : TO SAVE 3-D CLOUD FRACTION                     !
445     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
446     ! 
447     !SELECT OPTIONS IN GRRAD
448     !
449            NP3D=5         ! 3: ferrier's microphysics cloud scheme (only stratiform cloud) (set iflagliq>0 in radsw_param.f and radlw_param.f)
450                           ! 4: zhao/carr/sundqvist microphysics cloud (not in use for NMMB)
451                           ! NEW OPTIONS FOR NMMB
452                           ! 5: NAM stratiform + convective cloud optical depth and fraction (set iflagliq=0 in radsw_param.f and radlw_param.f)
453            ISOL=0         ! 0: use a fixed solar constant value (default)
454                           ! 1: use 11-year cycle solar constant table
455            ICO2=0         ! 0: use prescribed global mean co2   (default)
456                           ! 1: use observed co2 annual mean value only
457                           ! 2: use obs co2 monthly data with 2-d variation
458            ICWP=1         ! control flag for cloud generation schemes
459                           ! 0: use diagnostic cloud scheme
460                           ! 1: use prognostic cloud scheme (use with NMMB for NP3D=3,5)  
461            IALB=1         ! control flag for surface albedo schemes
462                           ! 0: climatology, based on surface veg types   !ONLY THIS ONE WORKS
463                           ! 1: modis retrieval based surface albedo scheme
464            IEMS=0         ! control flag for surface emissivity schemes 
465                           ! 0: fixed value of 1.0   (default)
466                           ! 1: varying value based on surface veg types
467            IAER= 11       ! flag for aerosols scheme selection (all options work for NMMB)
468                           ! - 3-digit aerosol flag (volc,lw,sw)     
469                           !    =  0: turn all aeros effects off (sw,lw,volc)                   !
470                           !    =  1: use clim tropspheric aerosol for sw only                  !
471                           !    = 10: use clim tropspheric aerosol for lw only                  !
472                           !    = 11: use clim tropspheric aerosol for both sw and lw           !
473                           !    =100: volc aerosol only for both sw and lw                      !
474                           !    =101: volc and clim trops aerosol for sw only                   !
475                           !    =110: volc and clim trops aerosol for lw only                   !
476                           !    =111: volc and clim trops aerosol for both sw and lw            !
477                           !    =  2: gocart/BSC-Dust tropspheric aerosol for sw only           !
478                           !    = 20: gocart/BSC-Dust tropspheric aerosol for lw only           !
479                           !    = 22: gocart/BSC-Dust tropspheric aerosol for both sw and lw    !
480                           !    =102: volc and gocart trops aerosol for sw only                 !
481                           !    =120: volc and gocart trops aerosol for lw only                 !
482                           !    =122: volc and gocart trops aerosol for both sw and lw          !
483            FHAER=0.       ! = 0 aerosol determined from gocart clim  !!Only this option works since on-line aerosol is not available in this version
484                           ! > 0. and <99999. aerosol determined from fcst and gocart clim
485                           ! >99999. aerosol determined from fcst
486     
487     
488     !
489            NFXR=27        ! second dimension of input/output array fluxr
490            NTRAC=3        ! dimension veriable for array oz
491            LSSAV=.FALSE.  ! logical flag for store 3-d cloud field
492            LDIAG3D=.FALSE.! logical flag for store 3-d diagnostic fields
493            LPRNT=.FALSE.
494            SASHAL=.FALSE. ! New Massflux based shallow convection           ! Not in use for NMMB
495            NTOZ=0         !  =0 climatological ozone profile                !
496                           !  >0 interactive ozone profile                   !
497            NCLDX=1        !  : only used when ntcw .gt. 0
498            NTCW=3         !  =0 no cloud condensate calculated              !
499                           !  >0 array index location for cloud condensate   !
500            IOVR_SW=1      ! 0 sw: random overlap clouds
501                           ! 1 sw: max-random overlap clouds
502            IOVR_LW=1      ! 0 lw: random overlap clouds
503                           ! 1 lw: max-random overlap clouds
504     !
505     !ADDITIONAL OPTIONS FOR NMMB
506     !
507            ALBTYPE=0      ! 0: CALCULATES ALBEDO FROM NMMB MONTHLY CLIMATOLOGY AS IN GFDL RADIATION
508     !
509            FHSWR=FLOAT(NRADS*DT_INT)/3600.   ! [h]
510            FHLWR=FLOAT(NRADL*DT_INT)/3600.   ! [h]
511            DTLW =FLOAT(NRADL*DT_INT)         ! [s]
512            DTSW =FLOAT(NRADS*DT_INT)         ! [s]
513            LSSWR=MOD(NTIMESTEP,NRADS)==0
514            LSLWR=MOD(NTIMESTEP,NRADL)==0
515     !
516     !INIT RADIATION
517     !
518            DO L=1,LM+1
519             RSGM(L)=SGM(L)
520            ENDDO
521     
522            CALL n_RADINIT (RSGM, LM, IFLIP, NP3D, ISOL, ICO2,     &
523                           ICWP, IALB, IEMS, IAER, JDAT, MYPE,   &
524                           FHSWR, FHLWR, FHAER )
525     !
526            LONR             =ITE-ITS+1
527            LATS_NODE_R      =JTE-JTS+1
528            LATR             =1                                 
529            IPT_LATS_NODE_R  =1
530     !
531            DO J=JTS,JTE
532            GLOBAL_LATS_R(J)= J-JTS+1
533            LONSPERLAR(J)   = ITE-ITS+1
534            SINLAT_R(J)     = SIN(GLAT( (ITS+ITE)/2 ,J))
535            COSLAT_R(J)     = SQRT( 1.d0 - SINLAT_R(J)*SINLAT_R(J) )
536            DO I=ITS,ITE
537            XLON(I,J)      = GLON(I,J)
538            SFCALBEDO(I,J) = ALBEDO(I,J)
539            ENDDO
540            ENDDO
541     !
542            CALL n_ASTRONOMY                                                   &
543     !  --- inputs:
544                ( LONSPERLAR, GLOBAL_LATS_R, SINLAT_R, COSLAT_R, XLON,       &
545                  FHSWR, JDAT, NRADS,                                        &
546                  LONR, LATS_NODE_R, LATR, IPT_LATS_NODE_R, LSSWR,           &
547     !  --- outputs:
548                  SOLCON, SLAG, SDEC, CDEC, COSZEN, COSZDG                   &
549                 )
550     ! ----
551     !
552     !OZONE CLIMATOLOGY
553     !from gfs_physics_initialize_mod.f
554     !
555             IF (NTOZ .LE. 0) THEN      ! DIAGNOSTIC OZONE, ONLY THIS ONE WORKS
556             LEVOZC  = 17
557             LATSOZC = 18
558             BLATC   = -85.0
559             TIMEOZC = 12
560             LATSOZP   = 2
561             LEVOZP    = 1
562             TIMEOZ    = 1
563             PL_COEFF  = 0
564             ENDIF
565             DPHIOZC = -(BLATC+BLATC)/(LATSOZC-1)
566     !
567     !from gloobr.f
568     !
569            IF (NTOZ .LE. 0) THEN      ! CLIMATOLOGICAL OZONE!
570     !
571              IDAY   = JDAT(3)
572              IMON   = JDAT(2)
573              MIDMON = DAYS(IMON)/2 + 1
574              CHANGE = FIRST .OR.( (IDAY .EQ. MIDMON) .AND. (JDAT(5).EQ.0) )
575     !
576              IF (CHANGE) THEN
577                  IF (IDAY .LT. MIDMON) THEN
578                      K1OZ = MOD(IMON+10,12) + 1
579                      MIDM = DAYS(K1OZ)/2 + 1
580                      K2OZ = IMON
581                      MIDP = DAYS(K1OZ) + MIDMON
582                   ELSE
583                      K1OZ = IMON
584                      MIDM = MIDMON
585                      K2OZ = MOD(IMON,12) + 1
586                      MIDP = DAYS(K2OZ)/2 + 1 + DAYS(K1OZ)
587                   ENDIF
588               ENDIF
589     !
590               IF (IDAY .LT. MIDMON) THEN
591                  ID = IDAY + DAYS(K1OZ)
592               ELSE
593                  ID = IDAY
594               ENDIF
595               FACOZ = REAL (ID-MIDM) / REAL (MIDP-MIDM)
596     !
597            ENDIF
598     !
599     !CLOUDS
600     !
601     !----------------------CONVECTION--------------------------------------
602     !  NRADPP IS THE NUMBER OF TIME STEPS TO ACCUMULATE CONVECTIVE PRECIP
603     !     FOR RADIATION
604     !   NOTE: THIS WILL NOT WORK IF NRADS AND NRADL ARE DIFFERENT UNLESS
605     !         THEY ARE INTEGER MULTIPLES OF EACH OTHER
606     !  CLSTP IS THE NUMBER OF HOURS OF THE ACCUMULATION PERIOD
607     !
608           NTSPH=NINT(3600./FLOAT(DT_INT))
609           NRADPP=MIN(NRADS,NRADL)
610           CLSTP=1.0*NRADPP/NTSPH
611           CONVPRATE=CUPRATE/CLSTP
612     !
613           LM1=LM-1
614     !
615           DO J=JTS,JTE
616           DO I=ITS,ITE
617     !
618             P8W(I,J,1)=PT
619     !
620             DO K=1,LM
621               P8W(I,J,K+1)=P8W(I,J,K)+PDSG1(K)+DSG2(K)*PD(I,J)
622               P_PHY(I,J,K)=SGML2(K)*PD(I,J)+PSGML1(K)
623               CCMID(I,J,K)=0.
624               CSMID(I,J,K)=0.
625             ENDDO
626           ENDDO
627           ENDDO
628     !
629           DO K=1,LM
630           DO J=JTS,JTE
631           DO I=ITS,ITE
632             CLDFRA(I,J,K)=0.
633             TAUTOTAL(I,J,K)=0.
634           ENDDO
635           ENDDO
636           ENDDO
637     !
638           DO J=JTS,JTE
639           DO I=ITS,ITE
640               CFRACH(I,J)=0.
641               CFRACL(I,J)=0.
642               CFRACM(I,J)=0.
643               CZMEAN(I,J)=0.
644               SIGT4(I,J)=0.
645           ENDDO
646           ENDDO
647     !
648           DO K=1,3
649           DO J=JTS,JTE
650           DO I=ITS,ITE
651             CLDCFR(I,J,K)=0.
652             MTOP(I,J,K)=0
653             MBOT(I,J,K)=0
654           ENDDO
655           ENDDO
656           ENDDO
657     !
658           DO J=JTS,JTE
659           DO I=ITS,ITE
660            CUTOP(I,J)=LM+1-HTOP(I,J)
661            CUBOT(I,J)=LM+1-HBOT(I,J)
662           ENDDO
663           ENDDO
664     !      
665     !-----------------------------------------------------------------------
666     !---  COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION  (Ferrier, Nov '04)
667     !
668     !--- Assumes Gaussian-distributed probability density functions (PDFs) for
669     !    total relative humidity (RHtot) within the grid for convective and
670     !    grid-scale cloud processes.  The standard deviation of RHtot is assumed
671     !    to be larger for convective clouds than grid-scale (stratiform) clouds.
672     !-----------------------------------------------------------------------
673     !
674            DO J=JTS,JTE
675            DO I=ITS,ITE 
676     !
677             DO 255 L=1,LM
678     !
679                 WV=MAX(EPSQ,Q(I,J,L))/(1.-MAX(EPSQ,Q(I,J,L)))               !--- Water vapor mixing ratio
680     !
681                 QICE=MAX(WATER(I,J,L,P_QS),0.)
682     !
683                 QCLD=QICE+MAX(WATER(I,J,L,P_QC),0.)                         !--- Total cloud water + ice mixing ratio
684     !
685                 IF (QCLD .LE. EPSQ) GO TO 255                               !--- Skip if no condensate is present
686                 CLFR=H0
687     !
688                 WV=MAX(EPSQ,Q(I,J,L))/(1.-MAX(EPSQ,Q(I,J,L)))
689     !
690     !--- Saturation vapor pressure w/r/t water ( >=0C ) or ice ( <0C )
691     !
692                 ESAT=1000.*FPVS(T(I,J,L))                                   !--- Saturation vapor pressure (Pa)
693                 QSAT=EP_2*ESAT/(P_PHY(I,J,L)-ESAT)                          !--- Saturation mixing ratio
694     !
695                 RHUM=WV/QSAT                                                !--- Relative humidity
696     !
697     !--- Revised cloud cover parameterization (temporarily ignore rain)
698     !
699                 RHtot=(WV+QCLD)/QSAT                                        !--- Total relative humidity
700     !
701                 LCNVT=NINT(CUTOP(I,J))
702                 LCNVT=MIN(LM,LCNVT)
703                 LCNVB=NINT(CUBOT(I,J))
704                 LCNVB=MIN(LM,LCNVB)
705                 IF (L.GE.LCNVT .AND. L.LE.LCNVB) THEN
706                    SDM=CVSDM
707                 ELSE
708                    SDM=STSDM
709                 ENDIF
710                 ARG=(RHtot-RHgrd)/SDM
711                 IF (ARG.LE.DXSD2 .AND. ARG.GE.DXSD2N) THEN
712                    CLFR=HALF
713                 ELSE IF (ARG .GT. DXSD2) THEN
714                    IF (ARG .GE. XSDmax) THEN
715                       CLFR=H1
716                    ELSE
717                       IXSD=INT(ARG/DXSD+HALF)
718                       IXSD=MIN(NXSD, MAX(IXSD,1))
719                       CLFR=HALF+AXSD(IXSD)
720                    ENDIF              !--- End IF (ARG .GE. XSDmax)
721                 ELSE
722                    IF (ARG .LE. XSDmin) THEN
723                       CLFR=H0
724                    ELSE
725                       IXSD=INT(ARG/DXSD1+HALF)
726                       IXSD=MIN(NXSD, MAX(IXSD,1))
727                       CLFR=HALF-AXSD(IXSD)
728                       IF (CLFR .LT. CLFRmin) CLFR=H0
729                    ENDIF        !--- End IF (ARG .LE. XSDmin)
730                 ENDIF           !--- IF (ARG.LE.DXSD2 .AND. ARG.GE.DXSD2N)
731                 CSMID(I,J,L)=CLFR
732     !
733     255       CONTINUE         !--- End DO L=1,LML
734     
735            ENDDO ! End DO I=ITS,ITE
736            ENDDO ! End DO J=JTS,JTE
737     
738     !***********************************************************************
739     !******************  END OF GRID-SCALE CLOUD FRACTIONS  ****************
740     !
741     !---  COMPUTE CONVECTIVE CLOUD COVER FOR RADIATION
742     !
743     !--- The parameterization of Slingo (1987, QJRMS, Table 1, p. 904) is
744     !    used for convective cloud fraction as a function of precipitation
745     !    rate.  Cloud fractions have been increased by 20% for each rainrate
746     !    interval so that shallow, nonprecipitating convection is ascribed a
747     !    constant cloud fraction of 0.1  (Ferrier, Feb '02).
748     !***********************************************************************
749     !
750           IF (CNCLD) THEN
751     
752            DO J=JTS,JTE
753             DO I=ITS,ITE
754     !
755     !***  CLOUD TOPS AND BOTTOMS COME FROM CUCNVC
756     !     Convective clouds need to be at least 2 model layers thick
757     !
758               IF (CUBOT(I,J)-CUTOP(I,J) .GT. 1.0) THEN
759     !--- Compute convective cloud fractions if appropriate  (Ferrier, Feb '02)
760                 CLFR=CC(1)
761                 PMOD=CUPPT(I,J)*CONVPRATE
762                 IF (PMOD .GT. PPT(1)) THEN
763                   DO NC=1,10
764                     IF(PMOD.GT.PPT(NC)) NMOD=NC
765                   ENDDO
766                   IF (NMOD .GE. 10) THEN
767                     CLFR=CC(10)
768                   ELSE
769                     CC1=CC(NMOD)
770                     CC2=CC(NMOD+1)
771                     P1=PPT(NMOD)
772                     P2=PPT(NMOD+1)
773                     CLFR=CC1+(CC2-CC1)*(PMOD-P1)/(P2-P1)
774                   ENDIF      !--- End IF (NMOD .GE. 10) ...
775                   CLFR=MIN(H1, CLFR)
776                 ENDIF        !--- End IF (PMOD .GT. PPT(1)) ...
777     !
778     !***  ADD LVL TO BE CONSISTENT WITH OTHER WORKING ARRAYS
779     !
780                 LCNVT=NINT(CUTOP(I,J))
781                 LCNVT=MIN(LM,LCNVT)
782                 LCNVB=NINT(CUBOT(I,J))
783                 LCNVB=MIN(LM,LCNVB)
784     !
785     !--- Build in small amounts of subgrid-scale convective condensate
786     !    (simple assumptions), but only if the convective cloud fraction
787     !    exceeds that of the grid-scale cloud fraction
788     !
789                 DO L=LCNVT,LCNVB
790                   ARG=MAX(H0, H1-CSMID(I,J,L))
791                   CCMID(I,J,L)=MIN(ARG,CLFR)
792                 ENDDO           !--- End DO LL=LCNVT,LCNVB
793               ENDIF             !--- IF (CUBOT(I,J)-CUTOP(I,J) .GT. 1.0) ...
794             ENDDO               ! End DO I=ITS,ITE
795            ENDDO                ! End DO J=JTS,JTE
796           ENDIF                 !--- End IF (CNCLD) ...
797     !
798     !*********************************************************************
799     !***************  END OF CONVECTIVE CLOUD FRACTIONS  *****************
800     !*********************************************************************
801     !***
802     !*** INITIALIZE ARRAYS FOR USES LATER
803     !***
804     
805           DO I=ITS,ITE
806           DO J=JTS,JTE
807     !
808           LML=LM
809     !***
810     !*** NOTE: LAYER=1 IS THE SURFACE, AND LAYER=2 IS THE FIRST CLOUD
811     !***       LAYER ABOVE THE SURFACE AND SO ON.
812     !***
813           KTOP(I,J,1)=LM+1
814           KBTM(I,J,1)=LM+1
815           CAMT(I,J,1)=1.0
816           KCLD(I,J)=2
817     !
818           DO 510 L=2,LM+1
819           CAMT(I,J,L)=0.0
820           KTOP(I,J,L)=1
821           KBTM(I,J,L)=1
822       510 CONTINUE
823     !### End changes so far
824     !***
825     !*** NOW CALCULATE THE AMOUNT, TOP, BOTTOM AND TYPE OF EACH CLOUD LAYER
826     !*** CLOUD TYPE=1: STRATIFORM CLOUD
827     !***       TYPE=2: CONVECTIVE CLOUD
828     !*** WHEN BOTH CONVECTIVE AND STRATIFORM CLOUDS EXIST AT THE SAME POINT,
829     !*** SELECT CONVECTIVE CLOUD WITH THE HIGHER CLOUD FRACTION.
830     !*** CLOUD LAYERS ARE SEPARATED BY TOTAL ABSENCE OF CLOUDINESS.
831     !*** NOTE: THERE IS ONLY ONE CONVECTIVE CLOUD LAYER IN ONE COLUMN.
832     !*** KTOP AND KBTM ARE THE TOP AND BOTTOM OF EACH CLOUD LAYER IN TERMS
833     !*** OF MODEL LEVEL.
834     !***
835           NEW_CLOUD=.TRUE.
836     !
837           DO L=2,LML
838             LL=LML-L+1                                  !-- Model layer
839             CLFR=MAX(CCMID(I,J,LL),CSMID(I,J,LL))       !-- Cloud fraction in layer
840             CLFR1=MAX(CCMID(I,J,LL+1),CSMID(I,J,LL+1))  !-- Cloud fraction in lower layer
841     !-------------------
842             IF (CLFR .GE. CLFRMIN) THEN
843     !--- Cloud present at level
844               IF (NEW_CLOUD) THEN
845     !--- New cloud layer
846                 IF(L==2.AND.CLFR1>=CLFRmin)THEN
847                   KBTM(I,J,KCLD(I,J))=LL+1
848                   CAMT(I,J,KCLD(I,J))=CLFR1
849                 ELSE
850                   KBTM(I,J,KCLD(I,J))=LL
851                   CAMT(I,J,KCLD(I,J))=CLFR
852                 ENDIF
853                 NEW_CLOUD=.FALSE.
854               ELSE
855     !--- Existing cloud layer
856                 CAMT(I,J,KCLD(I,J))=AMAX1(CAMT(I,J,KCLD(I,J)), CLFR)
857               ENDIF        ! End IF (NEW_CLOUD .EQ. 0) ...
858             ELSE IF (CLFR1 .GE. CLFRMIN) THEN
859     !--- Cloud is not present at level but did exist at lower level, then ...
860               IF (L .EQ. 2) THEN
861     !--- For the case of ground fog
862                KBTM(I,J,KCLD(I,J))=LL+1
863                CAMT(I,J,KCLD(I,J))=CLFR1
864               ENDIF
865               KTOP(I,J,KCLD(I,J))=LL+1
866               NEW_CLOUD=.TRUE.
867               KCLD(I,J)=KCLD(I,J)+1
868               CAMT(I,J,KCLD(I,J))=0.0
869             ENDIF
870     !-------------------
871           ENDDO      !--- End DO L loop
872     !***
873     !*** THE REAL NUMBER OF CLOUD LAYERS IS (THE FIRST IS THE GROUND;
874     !*** THE LAST IS THE SKY):
875     !***
876           NCLDS(I,J)=KCLD(I,J)-2
877           NCLD=NCLDS(I,J)
878     !***
879     !***  NOW CALCULATE CLOUD RADIATIVE PROPERTIES
880     !***
881           IF(NCLD.GE.1)THEN
882     !***
883     !*** NOTE: THE FOLLOWING CALCULATIONS, THE UNIT FOR PRESSURE IS MB!!!
884     !***
885             DO 580 NC=2,NCLD+1
886     !
887             TauC=0.    !--- Total optical depth for each cloud layer (solar & longwave)
888             QSUM=0.0
889             NKTP=LM+1
890             NBTM=0
891             BITX=CAMT(I,J,NC).GE.CLFRMIN
892             NKTP=MIN(NKTP,KTOP(I,J,NC))
893             NBTM=MAX(NBTM,KBTM(I,J,NC))
894     !
895             DO LL=NKTP,NBTM
896               L=NBTM-LL+NKTP 
897               IF(LL.GE.KTOP(I,J,NC).AND.LL.LE.KBTM(I,J,NC).AND.BITX)THEN
898                 PRS1=P8W(I,J,L)*0.01 
899                 PRS2=P8W(I,J,L+1)*0.01
900                 DELP=PRS2-PRS1
901                 TCLD=T(I,J,L)-T0C 
902                 QSUM=QSUM+Q(I,J,L)*DELP*(PRS1+PRS2)      & 
903          &           /(120.1612*SQRT(T(I,J,L)))
904     !
905                 CTau=0.
906     !-- For crude estimation of convective cloud optical depths
907                 IF (CCMID(I,J,L) .GE. CLFRmin) THEN
908                   IF (TCLD .GE. TRAD_ice) THEN
909                     CTau=CTauCW            !--- Convective cloud water
910                   ELSE
911                     CTau=CTauCI            !--- Convective ice
912                   ENDIF
913                 ENDIF
914     !
915     !-- For crude estimation of grid-scale cloud optical depths
916     !
917     !--   => The following 2 lines were intended to reduce cloud optical depths further
918     !        than what's parameterized in the NAM and what's theoretically justified
919                 CTau=CTau+ABSCOEF_W*WATER(I,J,L,P_QC)+ABSCOEF_I*WATER(I,J,L,P_QS)
920     
921                 TAUTOTAL(I,J,L)=CTau*DELP                          !Total model level cloud optical depth
922                 CLDFRA(I,J,L)=MAX(CCMID(I,J,LL),CSMID(I,J,LL))     !Cloud fraction at model level           
923                 TauC=TauC+DELP*CTau                                !Total cloud optical depth as in GFDL
924     !
925               ENDIF      !--- End IF(LL.GE.KTOP(I,NC) ....
926             ENDDO        !--- End DO LL
927     
928     !
929       580   CONTINUE
930     !
931           ENDIF        !!!NCLD.GE.1
932     !
933           ENDDO  !  DO I=ITS,ITE
934           ENDDO  !  DO J=JTS,JTE
935     !
936     ! Main domain loop: calling grrad
937     !
938           FLGMIN_L(1)= 0.20d0 ! --- for ferrier
939           FLUXR_V= 0.0d0      ! to save 2-d fields
940     
941           CV (1)=0.d0         ! not in use 
942           CVB(1)=0.d0         ! not in use
943           CVT(1)=0.d0         ! not in use
944     !
945           DO J=JTS,JTE  !start grrad loop column by column
946           DO I=ITS,ITE
947     !
948            CZMEAN(I,J)=COSZEN(I,J)
949            XLAT(1)=GLAT(I,J)
950            TSEA(1)=TSKIN(I,J)
951            TISFC(1)=TSKIN(I,J)                  ! change later if necessary
952            ZORL(1)=Z0(I,J)*100.d0
953            SNWDPH(1)=SI(I,J)                    ! snwdph[mm]
954            SNCOVR(1)=SNOW(I,J)/(SNOW(I,J)+70.)  ! FORMULATION OF MARSHALL ET AL. 1994
955            SNOALB(1)=MXSNAL(I,J)
956            HPRIME_V(1)=STDH(I,J)
957     
958            IF(SICE(I,J).GT.0.5) THEN              ! slmsk - ocean  - 0
959             SLMSK(1)= 2.0d0                       !         land   - 1
960             FICE(1)=SICE(I,J)                     ! change this later
961            ELSE                                   !         seaice - 2
962             SLMSK(1)= 1.0d0-SM(I,J)               !
963             FICE(1)= 0.0d0                        ! change this later
964            ENDIF
965     !
966     !!!ALBEDOS
967     ! 
968            IF (ALBTYPE==0) THEN
969     !..... THE FOLLOWING CODE GETS ALBEDO FROM PAYNE,1972 TABLES IF
970     !         1) OPEN SEA POINT (SLMSK=1);2) KALB=0
971                 IQ=INT(TWENTY*HP537+ONE)
972                 IF(CZMEAN(I,J).GT.0.0 .AND. SM(I,J).GT.0.5) THEN
973                     ZEN=DEGRAD1*ACOS(MAX(CZMEAN(I,J),0.0))
974                     IF(ZEN.GE.H74E1) JX=INT(HAF*(HNINETY-ZEN)+ONE)
975                     IF(ZEN.LT.H74E1.AND.ZEN.GE.FIFTY) &
976                     JX=INT(QUARTR*(H74E1-ZEN)+HNINE)
977                  IF(ZEN.LT.FIFTY) JX=INT(HP1*(FIFTY-ZEN)+H15E1)
978                  DZEN=-(ZEN-ZA(JX))/DZA(JX)
979                  ALB1=ALBD(IQ,JX)+DZEN*(ALBD(IQ,JX+1)-ALBD(IQ,JX))
980                  ALB2=ALBD(IQ+1,JX)+DZEN*(ALBD(IQ+1,JX+1)-ALBD(IQ+1,JX))
981                  SFCALBEDO(I,J)=ALB1+TWENTY*(ALB2-ALB1)*(HP537-TRN(IQ))
982                 ENDIF
983     !.....     VISIBLE AND NEAR IR DIFFUSE ALBEDO
984                  ALVD(1) = SFCALBEDO(I,J)
985                  ALND(1) = SFCALBEDO(I,J)
986     !.....     VISIBLE AND NEAR IR DIRECT BEAM ALBEDO
987                  ALVB(1) = SFCALBEDO(I,J)
988                  ALNB(1) = SFCALBEDO(I,J)
989     !
990     !--- Remove diurnal variation of land surface albedos (Ferrier, 6/28/05)
991     !--- Turn back on to mimic NAM 8/17/05
992     !
993     !.....     VISIBLE AND NEAR IR DIRECT BEAM ALBEDO, IF NOT OCEAN NOR SNOW
994     !        ..FUNCTION OF COSINE SOLAR ZENITH ANGLE..
995                 IF (SM(I,J).LT.0.5) THEN
996                  IF (SFCALBEDO(I,J).LE.0.5) THEN
997                  ALBD0=-18.0 * (0.5 - ACOS(CZMEAN(I,J))/PI)
998                  ALBD0=EXP (ALBD0)
999                  ALVD1=(ALVD(1) - 0.054313) / 0.945687
1000                  ALND1=(ALND(1) - 0.054313) / 0.945687
1001                  ALVB(1)=ALVD1 + (1.0 - ALVD1) * ALBD0
1002                  ALNB(1)=ALND1 + (1.0 - ALND1) * ALBD0
1003     ! !-- Put in an upper limit on beam albedos
1004                  ALVB(1)=MIN(0.5,ALVB(1))
1005                  ALNB(1)=MIN(0.5,ALNB(1))
1006                 END IF
1007                END IF
1008     
1009     !!! WE INTRODUCE HERE DIRECT AND DIFFUSE ALBEDO... FOR THIS OPTION, THERE IS A CHANGE IN GRRAD.f
1010                 ALVSF1(1)=ALVB(1) !For this option ALVSF1 is direct visible albedo
1011                 ALNSF1(1)=ALNB(1) !For this option ALNSF1 is direct nir albedo
1012                 ALVWF1(1)=ALVD(1) !For this option ALVWF1 is diffuse visible albedo
1013                 ALNWF1(1)=ALND(1) !For this option ALNWF1 is diffuse nir albedo
1014                 FACSF1(1)=0.      !not used with this option
1015                 FACWF1(1)=0.      !not used for this option
1016     !
1017            ENDIF
1018     !
1019     !---
1020           PRSI(1)=P8W(I,J,1)/1000.                                  ! [kPa]
1021     !
1022           DO L=1,LM
1023             PRSI(L+1)=P8W(I,J,L+1)/1000.                          ! (pressure on interface) [kPa]
1024             PRSL(L)=P_PHY(I,J,L)/1000.                            ! (pressure on mid-layer) [kPa] 
1025             PRSLK(L)=(PRSL(L)*0.01d0)**(R/CP)
1026             RTvR=1./(R*(Q(I,J,L)*0.608+1.-CW(I,J,L))*T(I,J,L))
1027             VVEL(L)=OMGALF(I,J,L)*1000.d0*PRSL(L)*RTvR            !not used
1028             GT(L)=T(I,J,L)
1029             GQ(L)=Q(I,J,L)
1030     !
1031             if(ntoz.le.0) then
1032               gr1(1,l,1)=0.d0
1033             else
1034               gr1(1,l,1)=max(o3(i,j,l),epso3)
1035             endif
1036     !
1037             GR1(1,L,2)=0.d0
1038             GR1(1,L,3)=CW(I,J,L)
1039             CLDCOV_V(L)=0.                                        !not used
1040             F_ICEC(L)=F_ICE(I,J,L)
1041             F_RAINC(L)=F_RAIN(I,J,L)
1042             R_RIME(L)=F_RIMEF(I,J,L)
1043             TAUCLOUDS(L)=TAUTOTAL(I,J,L)    !CLOUD OPTICAL DEPTH
1044             CLDF(L)=CLDFRA(I,J,L)           !CLOUD FRACTION 
1045           ENDDO
1046     !
1047     !---
1048           CALL n_GRRAD                                                             &
1049     !  ---  INPUTS:
1050                ( PRSI,PRSL,PRSLK,GT,GQ,GR1,VVEL,SLMSK,                           &
1051                  XLON(I,J),XLAT,TSEA,SNWDPH,SNCOVR,SNOALB,ZORL,HPRIME_V,         &
1052                  ALVSF1,ALNSF1,ALVWF1,ALNWF1,FACSF1,FACWF1,FICE,TISFC,           &
1053                  SOLCON,COSZEN(I,J),COSZDG(I,J),K1OZ,K2OZ,FACOZ,                 &
1054                  CV,CVT,CVB,IOVR_SW, IOVR_LW, F_ICEC, F_RAINC, R_RIME, FLGMIN_L, &
1055                  NP3D,NTCW,NCLDX,NTOZ,NTRAC,NFXR,                                &
1056                  DTLW,DTSW,LSSWR,LSLWR,LSSAV,LDIAG3D,SASHAL,                     &
1057                  1, 1, LM, IFLIP, MYPE, LPRNT,                                   &
1058     !  ---  ADDITIONAL INPUTS:
1059                  LATSOZC,LEVOZC,BLATC,DPHIOZC,TIMEOZC,                           &
1060                  TAUCLOUDS,CLDF,ALBTYPE,                                         &
1061     !  ---  OUTPUTS:
1062                  SWH,SFCNSW,SFCDSW,SFALB,                                        &
1063                  HLW,SFCDLW,TSFLW,                                               &
1064     !  ---  ADDITIONAL OUTPUT
1065                  TOAUSW,TOADSW,SFCCDSW,TOAULW,SFCUSW,                    & 
1066     !  ---  INPUT/OUTPUT:
1067                  FLUXR_V,CLDCOV_V                                        &
1068                )
1069     !
1070           DO L=1,LM
1071             RLWTT(I,J,L)=HLW(L)
1072             RSWTT(I,J,L)=SWH(L)
1073           ENDDO
1074     !
1075           RLWIN(I,J)=SFCDLW(1)
1076           RSWIN(I,J)=SFCDSW(1)
1077           RSWINC(I,J)=SFCCDSW(1) 
1078           RSWOUT(I,J)=RSWIN(I,J)*SFALB(1)
1079           RLWTOA(I,J)=TOAULW(1)
1080           RSWTOA(I,J)=TOAUSW(1)
1081     !
1082           ENDDO     ! --- END I LOOP for grrad
1083           ENDDO     ! --- END J LOOP for grrad
1084     !
1085     !*** --------------------------------------------------------------------------
1086     !***  DETERMINE THE FRACTIONAL CLOUD COVERAGE FOR HIGH, MID
1087     !***  AND LOW OF CLOUDS FROM THE CLOUD COVERAGE AT EACH LEVEL
1088     !***
1089     !***  NOTE: THIS IS FOR DIAGNOSTICS ONLY!!!
1090     !***
1091     !***
1092     !
1093            DO J=JTS,JTE
1094            DO I=ITS,ITE
1095     !!
1096            DO L=0,LM
1097              CLDAMT(L)=0.
1098            ENDDO
1099     !!
1100     !!***  NOW GOES LOW, MIDDLE, HIGH
1101     !!
1102            DO 480 NLVL=1,3
1103            CLDMAX=0.
1104            MALVL=LM
1105            LLTOP=LM+1-LTOP(NLVL)   !!!!COMES FROM GFDL INIT
1106     !!***
1107     !!***  GO TO THE NEXT CLOUD LAYER IF THE TOP OF THE CLOUD-TYPE IN
1108     !!***  QUESTION IS BELOW GROUND OR IS IN THE LOWEST LAYER ABOVE GROUND.
1109     !!***
1110            IF(LLTOP.GE.LM)GO TO 480
1111     !!
1112            IF(NLVL.GT.1)THEN
1113              LLBOT=LM+1-LTOP(NLVL-1)-1
1114              LLBOT=MIN(LLBOT,LM1)
1115            ELSE
1116              LLBOT=LM1
1117            ENDIF
1118     !!
1119            DO 435 L=LLTOP,LLBOT
1120            CLDAMT(L)=AMAX1(CSMID(I,J,L),CCMID(I,J,L))
1121            IF(CLDAMT(L).GT.CLDMAX)THEN
1122              MALVL=L
1123              CLDMAX=CLDAMT(L)
1124            ENDIF
1125        435 CONTINUE
1126     !!*********************************************************************
1127     !! NOW, CALCULATE THE TOTAL CLOUD FRACTION IN THIS PRESSURE DOMAIN
1128     !! USING THE METHOD DEVELOPED BY Y.H., K.A.C. AND A.K. (NOV., 1992).
1129     !! IN THIS METHOD, IT IS ASSUMED THAT SEPERATED CLOUD LAYERS ARE
1130     !! RADOMLY OVERLAPPED AND ADJACENT CLOUD LAYERS ARE MAXIMUM OVERLAPPED.
1131     !! VERTICAL LOCATION OF EACH TYPE OF CLOUD IS DETERMINED BY THE THICKEST
1132     !! CONTINUING CLOUD LAYERS IN THE DOMAIN.
1133     !!*********************************************************************
1134            CL1=0.0
1135            CL2=0.0
1136            KBT1=LLBOT
1137            KBT2=LLBOT
1138            KTH1=0
1139            KTH2=0
1140     !!
1141            DO 450 LL=LLTOP,LLBOT
1142            L=LLBOT-LL+LLTOP
1143            BIT1=.FALSE.
1144            CR1=CLDAMT(L)
1145            BITX=(P8W(I,J,L).GE.PTOPC(NLVL+1)).AND.                           &
1146           &     (P8W(I,J,L).LT.PTOPC(NLVL)).AND.                             &
1147           &     (CLDAMT(L).GT.0.0)
1148            BIT1=BIT1.OR.BITX
1149            IF(.NOT.BIT1)GO TO 450
1150     !!***
1151     !!***  BITY=T: FIRST CLOUD LAYER; BITZ=T:CONSECUTIVE CLOUD LAYER
1152     !!***  NOTE:  WE ASSUME THAT THE THICKNESS OF EACH CLOUD LAYER IN THE
1153     !!***         DOMAIN IS LESS THAN 200 MB TO AVOID TOO MUCH COOLING OR
1154     !!***         HEATING. SO WE SET CTHK(NLVL)=200*E2. BUT THIS LIMIT MAY
1155     !!***         WORK WELL FOR CONVECTIVE CLOUDS. MODIFICATION MAY BE
1156     !!***         NEEDED IN THE FUTURE.
1157     !!***
1158            BITY=BITX.AND.(KTH2.LE.0)
1159            BITZ=BITX.AND.(KTH2.GT.0)
1160     !!
1161            IF(BITY)THEN
1162              KBT2=L
1163              KTH2=1
1164            ENDIF
1165     !!
1166            IF(BITZ)THEN
1167              KTOP1=KBT2-KTH2+1
1168              DPCL=P_PHY(I,J,KBT2)-P_PHY(I,J,KTOP1)
1169              IF(DPCL.LT.CTHK(NLVL))THEN
1170                KTH2=KTH2+1
1171              ELSE
1172                KBT2=KBT2-1
1173              ENDIF
1174            ENDIF
1175            IF(BITX)CL2=AMAX1(CL2,CR1)
1176     !!***
1177     !!*** AT THE DOMAIN BOUNDARY OR SEPARATED CLD LAYERS, RANDOM OVERLAP.
1178     !!*** CHOOSE THE THICKEST OR THE LARGEST FRACTION AMT AS THE CLD
1179     !!*** LAYER IN THAT DOMAIN.
1180     !!***
1181            BIT2=.FALSE.
1182            BITY=BITX.AND.(CLDAMT(L-1).LE.0.0.OR. &
1183                 P8W(I,J,L-1).LT.PTOPC(NLVL+1))
1184            BITZ=BITY.AND.CL1.GT.0.0
1185            BITW=BITY.AND.CL1.LE.0.0
1186            BIT2=BIT2.OR.BITY
1187            IF(.NOT.BIT2)GO TO 450
1188     !!
1189     !!
1190            IF(BITZ)THEN
1191              KBT1=INT((CL1*KBT1+CL2*KBT2)/(CL1+CL2))
1192              KTH1=INT((CL1*KTH1+CL2*KTH2)/(CL1+CL2))+1
1193              CL1=CL1+CL2-CL1*CL2
1194            ENDIF
1195     !!
1196            IF(BITW)THEN
1197              KBT1=KBT2
1198              KTH1=KTH2
1199              CL1=CL2
1200            ENDIF
1201     !!
1202            IF(BITY)THEN
1203              KBT2=LLBOT
1204              KTH2=0
1205              CL2=0.0
1206            ENDIF
1207       450 CONTINUE
1208     !
1209             CLDCFR(I,J,NLVL)=AMIN1(1.0,CL1)
1210             MTOP(I,J,NLVL)=MIN(KBT1,KBT1-KTH1+1)
1211             MBOT(I,J,NLVL)=KBT1
1212     
1213       480 CONTINUE
1214     
1215           ENDDO ! End DO I=ITS,ITE
1216           ENDDO ! End DO J=ITS,JTE
1217     !!
1218           DO J=JTS,JTE
1219           DO I=ITS,ITE
1220     
1221             CFRACL(I,J)=CLDCFR(I,J,1)
1222             CFRACM(I,J)=CLDCFR(I,J,2)
1223             CFRACH(I,J)=CLDCFR(I,J,3)
1224             IF(CNCLD)THEN
1225               CFSmax=0.   !-- Maximum cloud fraction (stratiform component)
1226               CFCmax=0.   !-- Maximum cloud fraction (convective component)
1227               DO L=1,LM
1228                 CFSmax=MAX(CFSmax, CSMID(I,J,L) )
1229                 CFCmax=MAX(CFCmax, CCMID(I,J,L) )
1230               ENDDO
1231               ACFRST(I,J)=ACFRST(I,J)+CFSmax
1232               NCFRST(I,J)=NCFRST(I,J)+1
1233               ACFRCV(I,J)=ACFRCV(I,J)+CFCmax
1234               NCFRCV(I,J)=NCFRCV(I,J)+1
1235             ELSE
1236       !--- Count only locations with grid-scale cloudiness, ignore convective clouds
1237       !    (option not used, but if so set to the total cloud fraction)
1238               CFRAVG=1.-(1.-CFRACL(I,J))*(1.-CFRACM(I,J))*(1.-CFRACH(I,J))
1239               ACFRST(I,J)=ACFRST(I,J)+CFRAVG
1240               NCFRST(I,J)=NCFRST(I,J)+1
1241             ENDIF
1242     
1243           ENDDO  !  DO I=ITS,ITE
1244           ENDDO  !  DO J=JTS,JTE
1245     !
1246     !-----------------------------------------------------------------------
1247     !***  LONGWAVE
1248     !-----------------------------------------------------------------------
1249     !
1250           IF(MOD(NTIMESTEP,NRADL)==0)THEN
1251     !.......................................................................
1252     !$omp parallel do                                                       &
1253     !$omp& private(i,j,k,kflip,tdum)
1254     !.......................................................................
1255             DO J=JTS,JTE
1256               DO I=ITS,ITE
1257     !
1258                 TDUM=T(I,J,LM)
1259                 SIGT4(I,J)=STBOLT*TDUM*TDUM*TDUM*TDUM
1260     !
1261               ENDDO
1262             ENDDO
1263     !.......................................................................
1264     !$omp end parallel do
1265     !.......................................................................
1266     !
1267           ENDIF
1268     !      
1269           IF (FIRST) FIRST=.FALSE.
1270     
1271     !-----------------------------------------------------------------------
1272     !
1273           END SUBROUTINE RRTM
1274     !
1275     !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
1276     !-----------------------------------------------------------------------
1277           SUBROUTINE RRTM_INIT(EMISS,SFULL,SHALF,PPTOP,                     &
1278          &                     JULYR,MONTH,IDAY,GMT,                        &
1279          &                     CO2TF,                                       &
1280          &                     IDS, IDE, JDS, JDE, KDS, KDE,                &
1281          &                     IMS, IME, JMS, JME, KMS, KME,                &
1282          &                     ITS, ITE, JTS, JTE, KTS, KTE              )
1283     !-----------------------------------------------------------------------
1284           IMPLICIT NONE
1285     !-----------------------------------------------------------------------
1286           INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                     &
1287          &                     ,IMS,IME,JMS,JME,KMS,KME                     &
1288          &                     ,ITS,ITE,JTS,JTE,KTS,KTE
1289           INTEGER,INTENT(IN) :: JULYR,MONTH,IDAY,CO2TF
1290           REAL,INTENT(IN) :: GMT,PPTOP
1291           REAL,DIMENSION(KMS:KME),INTENT(IN) :: SFULL, SHALF
1292           REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: EMISS
1293     !
1294           INTEGER :: I,IHRST,J,N
1295           REAL :: PCLD,XSD,PI,SQR2PI
1296           REAL :: SSLP=1013.25
1297           REAL, PARAMETER :: PTOP_HI=150.,PTOP_MID=350.,PTOP_LO=642.,       &
1298          &                   PLBTM=105000.
1299     !-----------------------------------------------------------------------
1300     !***********************************************************************
1301     !-----------------------------------------------------------------------
1302     !
1303     !***  INITIALIZE DIAGNOSTIC LOW,MIDDLE,HIGH CLOUD LAYER PRESSURE LIMITS.
1304     !
1305           LTOP(1)=0
1306           LTOP(2)=0
1307           LTOP(3)=0
1308     !
1309           DO N=1,KTE
1310             PCLD=(SSLP-PPTOP*10.)*SHALF(N)+PPTOP*10.
1311             IF(PCLD>=PTOP_LO)LTOP(1)=N
1312             IF(PCLD>=PTOP_MID)LTOP(2)=N
1313             IF(PCLD>=PTOP_HI)LTOP(3)=N
1314     !       PRINT *,N,PCLD,SHALF(N),PSTAR,PPTOP
1315           ENDDO
1316     !***
1317     !***  ASSIGN THE PRESSURES FOR CLOUD DOMAIN BOUNDARIES
1318     !***
1319           PTOPC(1)=PLBTM
1320           PTOPC(2)=PTOP_LO*100.
1321           PTOPC(3)=PTOP_MID*100.
1322           PTOPC(4)=PTOP_HI*100.
1323     !
1324     !***  FOR NOW, GFDL RADIATION ASSUMES EMISSIVITY = 1.0
1325     !
1326           DO J=JTS,JTE
1327           DO I=ITS,ITE
1328             EMISS(I,J) = 1.0
1329           ENDDO
1330           ENDDO
1331     !
1332     !---  Calculate the area under the Gaussian curve at the start of the
1333     !---  model run and build the look up table AXSD
1334     !
1335           PI=ACOS(-1.)
1336           SQR2PI=SQRT(2.*PI)
1337           RSQR=1./SQR2PI
1338           DO I=1,NXSD
1339             XSD=REAL(I)*DXSD
1340             AXSD(I)=GAUSIN(XSD)
1341           ENDDO
1342     !
1343     !-----------------------------------------------------------------------
1344           END SUBROUTINE RRTM_INIT
1345     !-----------------------------------------------------------------------
1346     !----------------------------------------------------------------------
1347     !
1348           REAL FUNCTION GAUSIN(xsd)
1349           REAL, PARAMETER :: crit=1.e-3
1350           REAL A1,A2,RN,B1,B2,B3,SUM,xsd
1351     !
1352     !  This function calculate area under the Gaussian curve between mean
1353     !  and xsd # of standard deviation (03/22/2004  Hsin-mu Lin)
1354     !
1355           a1=xsd*RSQR
1356           a2=exp(-0.5*xsd**2)
1357           rn=1.
1358           b1=1.
1359           b2=1.
1360           b3=1.
1361           sum=1.
1362           do while (b2 .gt. crit)
1363              rn=rn+1.
1364              b2=xsd**2/(2.*rn-1.)
1365              b3=b1*b2
1366              sum=sum+b3
1367              b1=b3
1368           enddo
1369           GAUSIN=a1*a2*sum
1370           RETURN
1371           END FUNCTION GAUSIN
1372     !
1373     !-----------------------------------------------------------------------
1374     !-----------------------------------------------------------------------
1375     !
1376           END MODULE MODULE_RA_RRTM
1377     !
1378     !-----------------------------------------------------------------------
1379