/[MITgcm]/MITgcm_contrib/jscott/igsm/src/radia.F
ViewVC logotype

Annotation of /MITgcm_contrib/jscott/igsm/src/radia.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (hide annotations) (download)
Mon Apr 23 21:20:18 2007 UTC (18 years, 3 months ago) by jscott
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +129 -63 lines
bring igsm atmos code up to date

1 jscott 1.1 #include "ctrparam.h"
2     SUBROUTINE RADIA 5001.
3     C**** 5002.
4     C**** THIS SUBROUTINES ADDS THE RADIATION HEATING TO THE TEMPERATURES 5003.
5     C**** 5004.
6    
7     #include "BD2G04.COM"
8     #include "chem_para"
9     #include "chem_com"
10    
11     COMMON U,V,T,P,Q 5006.
12     COMMON/ADDALB/BVSURFA,XVSURFA,BNSURFA,XNSURFA
13     COMMON/WORK1/CONV(IM0,JM0,LM0),PK(IM0,JM0,LM0),PREC(IM0,JM0),
14     & TPREC(IM0,JM0), 5007.
15     * COSZ1(IO0,JM0),COSZ2(IO0,JM0),COSZA(IO0,JM0), 5008.
16     * TRINCG(IO0,JM0),BTMPW(IO0,JM0),SNFS(IO0,JM0,4),TNFS(IO0,JM0,4), 5009.
17     * TRHRS(IO0,JM0,3),SRHRS(IO0,JM0,3),ALB(IO0,JM0,9) 5010.
18     COMMON/WORK2/CLDSS(IM0,JM0,LM0),CLDMC(IM0,JM0,LM0), 5011.
19     * TOTCLD(36) 5012.
20     DIMENSION TRNFP0(JM0),TRNFP1(JM0),ALBJ(JM0,9)
21     real ODATA2(JM0,2),GDATA2(JM0,14),BDATA2(JM0,2),FDATA2(JM0,2),
22     * RQT2(JM0,3)
23     common/SURRAD/TRSURF(JM0,4),SRSURF(JM0,4)
24     common/FORAERSOL/FORSULF,FORBC
25     logical FORSULF,FORBC
26     C COMMON/WORK4/ IS BEING USED BY THE RADIATION ROUTINES 5013.
27     C 5014.
28     C RADCOM: CONTROL/INPUT PARAMETERS 5015.
29     C 5016.
30     COMMON/RADCOM/VADATA(11,4,3),DGLAT(46),DGLON(72),TMINSR,FULGAS(18)5017.
31     A ,FRACSL,RATQSL,FOGTSL,PTLISO,TLGRAD,TKCICE,FGOLDU(18)5018.
32     B ,FLONO3,FRAYLE,FCLDTR,FCLDSR,FALGAE,FMARCL,FEMTRA(6) 5019.
33     C ,WETTRA,WETSRA,DMOICE,DMLICE,LICETK,NTRCE,FZASRA(6) 5020.
34     D ,ID5(5),ITR(4),IMG(2),ILG(2),LAPGAS,KWVCON,NORMS0,NV 5021.
35     E ,KEEPRH,KEEPAL,ISOSCT,IHGSCT,KFRACC,KGASSR,KAERSR 5022.
36     F ,MARCLD,LAYTOP,LMR,LMRP,JMLAT,IMLON,KFORCE,LASTVC 5023.
37     C 5024.
38     C BASIC RADCOM INPUT DATA 5025.
39     C 5026.
40     G ,PLE(40),HLB(40),TLB(40),TLT(40),TL(40),U0GAS(40,9) 5027.
41     H ,ULGAS(40,9),TRACER(40,4),RTAU(40),QL(40),RHL(40) 5028.
42     I ,POCEAN,PEARTH,POICE,PLICE,AGESN,SNOWE,SNOWOI,SNOWLI 5029.
43     J ,TGO,TGE,TGOI,TGLI,TS,WS,WEARTH,ZOICE,FSPARE(200) 5030.
44     K ,S0,COSZ,PVT(11),BXA(153),SRBXAL(15,2),FRC(5),LUXGAS 5031.
45     L ,JYEARR,JDAYR,JLAT,ILON,MEANAL,KALVIS,ISPARE(25),SGPS5032.
46     C 5033.
47     C BASIC RADCOM OUTPUT DATA 5034.
48     C 5035.
49     M ,TRDFLB(40),TRUFLB(40),TRNFLB(40),TRFCRL(40),TRSLCR 5036.
50     N ,SRDFLB(40),SRUFLB(40),SRNFLB(40),SRFHRL(40),SRSLHR 5037.
51     O ,SRIVIS,SROVIS,PLAVIS,SRINIR,SRONIR,PLANIR,SRXATM(4) 5038.
52     P ,SRDVIS,SRUVIS,ALBVIS,SRDNIR,SRUNIR,ALBNIR,FSRNFG(4) 5039.
53     Q ,SRTVIS,SRRVIS,SRAVIS,SRTNIR,SRRNIR,SRANIR,FTRUFG(4) 5040.
54     R ,TRDFGW,TRUFGW,TRUFTW,BTEMPW,TRDFSL,TRUFSL,DTRUFG(4) 5041.
55     S ,TRSLTS,TRSLTG,TRSLWV,TRSLBS,TTRUFG,LBOTCL,LTOPCL 5042.
56     DIMENSION COE(39) 5043.
57     LOGICAL POLE,DC25,HPRNT,WRCLD,CLDFEED 5044.
58     #if ( defined OCEAN_3D )
59 jscott 1.2 #include "AGRID.h"
60 jscott 1.1 #endif
61 jscott 1.3 dimension SWNET(jm0,2),SWIN(jm0,2)
62 jscott 1.1
63     #if ( defined CLM )
64 jscott 1.3 #include "CLM.h"
65 jscott 1.1 #endif
66     c
67    
68     common /ATCO2/atm_co2(jm0)
69     common/conprn/HPRNT
70     common/COMCLD/READGHG,PCLOUD,WRCLD,NWRCLD,NWRCL,INYEAR,JNDAY
71     &,CFAEROSOL,ALFA,CFBC,cfvolaer
72     common/ghstfor/GHSFALB,GHSF,ALBCF,FVOLADD,STRARFOR,S0FOR,CO2FOR,
73     & CO2F,ghostfv(LM0+1),ghostf(LM0+1,JM0),ghflux(LM0+1,JM0)
74     logical GHSFALB,GHSF,GSOEQ,STRARFOR,S0FOR,CO2FOR
75     common/eqgso/GSOEQ
76     common/cldfdb/coefcl(3),CLDFEED
77     common/aexpc/AEXP,ISTRT1
78     common/ SNOWALB/FRSNALB,III
79     common/ S0XR/S0RATE,CFS0X
80     common/ BACKGRGHG/GHGBGR(5)
81     COMMON/CO2TRND/ALFFOR,CO2TR,YEARGT,CO2IN,INYRAD
82     dimension CLDSSF(JM0,LM0),CLDMCF(JM0,LM0)
83     &,BSO4LAND(JM0),BSO4OCEAN(JM0),BSO4TOTAL(JM0)
84     &,bcodmn(JM0,LM0,12)
85     dimension DSWSRF(jm0),DLWSRF(jm0),DSWVIS(jm0),DSWNIR(jm0)
86 jscott 1.3 & ,ULWSRF(jm0)
87 jscott 1.1 integer PCLOUD
88 jscott 1.3 ! common/TSUR/TSURFC(JM0,0:13),TSURFT(JM0),TSURFD(JM0),DTSURF(JM0)
89     ! *,cfcld(JM0,3)
90     #include "TSRF.COM"
91 jscott 1.1 CHARACTER*4 JMNTHF,JMLAST
92     DATA JMLAST /'LAST'/
93     DATA TF/273.16/,TCIR/258.16/,STBO/.567257E-7/,IFIRST/1/,JDLAST/-9/5045.
94     DATA IRFIRST /1/
95     C **** CLEAR SKY
96     dimension SRHRCL(JM0),TRHRCL(JM0),ALBCL(JM0),SNP1CL(JM0),
97     *SNP0CL(JM0),TRINCL(JM0),TRP0CL(JM0),TRP1CL(JM0)
98     dimension RTAU0(LM0),TOTCLD0(LM0)
99     common/clrsk/CLEAR(JM0),NCLR(JM0),AJCLR(JM0,12),BJCLR(JM0,12),
100     * CJCLR(JM0,12)
101     integer CLEAR
102     C AJCLR
103     C 1 SW INC AT P0 RD (AJ(1))
104     C 2 SW ABS BELOW P0 RD (AJ(2))
105     C 3 SW ABS BELOW P1 RD (AJ(3))
106     C 4 SW ABS AT Z0 RD (AJ(6))
107     C 5 SW INC AT Z0 RD (AJ(5))
108     C 6 LW INC AT Z0 RD (AJ(67))
109     C 7 NET LW AT Z0 SF (AJ(9))
110     C 8 NET LW AT P0 RD (AJ(7))
111     C 9 NET LW AT P1 RD (AJ(8))
112     C 10 NET RAD AT P0 DG (AJ(10))
113     C 11 NET RAD AT P1 DG (AJ(11))
114     C 12 NET RAD AT Z0 DG (AJ(12))
115     C **** CLEAR SKY
116     C**** 5046.
117     C**** FDATA 2 LAND COVERAGE (1) 5047.
118     C**** 3 RATIO OF LAND ICE COVERAGE TO LAND COVERAGE (1) 5048.
119     C**** 5049.
120     C**** ODATA 1 OCEAN TEMPERATURE (C) 5050.
121     C**** 2 RATIO OF OCEAN ICE COVERAGE TO WATER COVERAGE (1) 5051.
122     C**** 5052.
123     C**** GDATA 1 OCEAN ICE SNOW AMOUNT (KG/M**2) 5053.
124     C**** 2 EARTH SNOW AMOUNT (KG/M**2) 5054.
125     C**** 3 OCEAN ICE TEMPERATURE OF FIRST LAYER (C) 5055.
126     C**** 4 EARTH TEMPERATURE OF FIRST LAYER (C) 5056.
127     C**** 5 EARTH WATER OF FIRST LAYER (KG/M**2) 5057.
128     C**** 6 EARTH ICE OF FIRST LAYER (KG/M**2) 5058.
129     C**** 11 AGE OF SNOW (DAYS) 5059.
130     C**** 12 LAND ICE SNOW AMOUNT (KG/M**2) 5060.
131     C**** 13 LAND ICE TEMPERATURE OF FIRST LAYER (C) 5061.
132     C**** 5062.
133     C**** BLDATA 1 COMPOSITE SURFACE WIND MAGNITUDE (M/S) 5063.
134     C**** 2 COMPOSITE SURFACE AIR TEMPERATURE (K) 5064.
135     C**** 5 FREE 5065.
136     C**** 5066.
137     C**** VDATA 1-8 EARTH RATIOS FOR THE 8 VEGETATION TYPES (1) 5067.
138     C**** 9 WATER FIELD CAPACITY OF FIRST LAYER (KG/M**2) 5068.
139     C**** 5069.
140     dimension GHGBGR0(5)
141    
142     IF(MODRD.EQ.0) IDACC(2)=IDACC(2)+1 5070.
143     IF (IFIRST.NE.1) GO TO 50 5071.
144     if(GHSF)then
145     do j=1,jm
146     ghflux(1,j)=ghostf(1,j)
147     do l=2,lm0+1
148     ghflux(l,j)= ghflux(l-1,j)+ghostf(l,j)
149     enddo
150     enddo
151     do l=1,lm0+1
152     print *,l,ghflux(l,1)
153     enddo
154     endif
155     JDAYR=JNDAY
156     JYEARR=INYEAR
157     c GHGBGR(1)=GHGBGR0(1)
158     c GHGBGR(2)=GHGBGR0(2)
159     c GHGBGR(3)=GHGBGR0(3)
160     c GHGBGR(4)=GHGBGR0(4)*1000.
161     c GHGBGR(5)=GHGBGR0(5)*1000.
162     c print *,'Background GHGs'
163     print '(5E12.4)',GHGBGR
164     nreadcld=0
165     nrbyyr=24*365/5
166     nrcldmax=20*nrbyyr
167     c print *,' CLOUDS for ',nrcldmax/nrbyyr,' years'
168     print *,'RADIA'
169     KTREND=-CO2
170     JDAY00=-1
171     if(CLDFEED)then
172     c print *,' radia.f coefcl=',coefcl
173     print *,' for low and middle clouds',coefcl(1)
174     print *,' for top clouds',coefcl(2)
175     print *,' for MC clouds',coefcl(3)
176     endif
177     print *,' READGHG=',READGHG
178     print *,' CFAEROSOL=',CFAEROSOL
179     #if ( defined PREDICTED_BC)
180 jscott 1.3 if(READGHG.eq.0)then
181     ! For runs with black carbon predicted by chemistry model
182 jscott 1.1 print *,'With black carbon forcing'
183 jscott 1.3 print *,' CFBC=',CFBC
184     nrdbc=0
185     read(769),bcodmn
186     endif
187 jscott 1.1 #endif
188     do j=1,jm
189     BSO4TOTAL(j)=1.0
190     BSO4LAND(j)=1.0
191     BSO4OCEAN(j)=1.0
192     enddo
193     print *,' separate caclulations for land and ocean'
194     DC25=.TRUE.
195     ! DC25=.FALSE.
196     if(DC25)then
197     print *,' with DC'
198     else
199     print *,' without DC'
200     print *,' subroutine COSZR'
201     end if
202     if(abs(PCLOUD-3.).gt.1.5.and..NOT.WRCLD)IFIRST=0 5072.
203     LMP1=LM+1 5072.1
204     DTCNDS=NCNDS*DT 5073.
205     C**** SET THE CONTROL PARAMETERS FOR THE RADIATION 5074.
206     JMLAT=JM 5074.1
207     if(JM.ne.24) then
208     DO J=1,JMLAT
209     DGLAT(J)=acos(COSP(J))*360./TWOPI
210     if(J.le.JMLAT/2)DGLAT(J)=-DGLAT(J)
211     END DO
212     endif
213     IMLON=IO 5074.2
214     LMR=LM+3 5075.
215     COEX=.01*GRAV*KAPA/RGAS 5076.
216     PSFMPT=PSF-PTOP 5077.
217     DO 30 L=1,LM 5078.
218     COE(L)=DTCNDS*COEX/DSIG(L) 5079.
219     30 PLE(L)=SIGE(L)*(PSF-PTOP)+PTOP 5080.
220     PLE(LMP1)=PTOP 5081.
221     PLE(LM+2)=.5*PTOP 5082.
222     PLE(LMR)=.2*PTOP 5083.
223     PLE(LMR+1)=1.E-5 5084.
224     DO 40 LR=LMP1,LMR 5085.
225     COE(LR)=DT*NRAD*COEX/(PLE(LR)-PLE(LR+1)) 5086.
226     QL(LR)=.3E-5 5087.
227     40 RTAU(LR)=0. 5088.
228     DPMICE=10. 5089.
229     C S0X=1. 5089.1
230     KSTREND=S0RATE*100.
231     print *,'S0RATE=',S0RATE,' KSTREND=',KSTREND
232     print *,'From first call of radia CFS0X=',CFS0X
233     #if ( defined FIXED_FOR )
234     print *,'Run with fixed forcing'
235     #endif
236     C
237     if(GSOEQ)then
238     if(KSTREND.ne.0)then
239     print *,'Wrong setting of GSOEQ and KSTREND'
240     print *,' GSOEQ=',GSOEQ,' KSTREND=',KSTREND
241     c stop
242     endif
243     if(YEARGT.ne.1860)then
244     print *,'Wrong setting of GSOEQ and YEARGT'
245     print *,' GSOEQ=',GSOEQ,' YEARGT=',YEARGT
246     stop
247     endif
248     C For an equilibrium run for 1860 conditions for GSO runs
249     print *,
250     & 'An equilibrium run for 1860 conditions for GSO runs'
251     S0X=1365.3596/1367.
252     print *,'Sconstant=',S0X*1367.
253     RVOL=0.012
254     FVOL=0.0045
255     FGOLDU(1)=(RVOL+FVOL)/RVOL
256     print *,'Additional STAER=',FVOL,
257     & ' FGOLDU(1)=',FGOLDU(1)
258     C
259     endif
260     if(CO2FOR)then
261     print *,' Forcing due to ',CO2F,'xCO2 is calculated',
262     * 'in stead of cloud forcing'
263     endif
264     if(S0FOR)then
265     print *,' Forcing due to change in solar constant ',
266     * 'is calculated in stead of cloud forcing'
267     endif
268     if(FORBC)then
269     print *,' Forcing due to black carbon ',
270     * 'is calculated in stead of cloud forcing'
271     endif
272     if(STRARFOR)then
273     RVOL=0.012
274     FVOL=FVOLADD
275     FGOLDU1=(RVOL+FVOL)/RVOL
276     print *,'Additional STAER=',FVOL,
277     & ' FGOLDU1=',FGOLDU1
278     endif
279     if(GHSFALB)then
280     Print *,' Run with ghost albedo forcing'
281     print *,'ALBCF=',ALBCF
282     endif
283     if(GHSF)then
284     Print *,' Run with ghost forcing'
285     do l=1,lm+1
286     ghostfv(l)=0.0
287     do j=1,jm
288     ghostfv(l)=ghostfv(l)+ghostf(l,j)
289     enddo
290     ghostfv(l)=ghostfv(l)/float(jm)
291     enddo
292     print *,' Ghost forcing DFS=',ghostfv(1)
293     do l=LM+1,2,-1
294     print *,' Forsing at laeyr',l-1,'=',ghostfv(l)
295     enddo
296     c print *,' 2D ghost forcing'
297     endif
298     S0X=CFS0X*S0X
299     S0X0=S0X
300 jscott 1.3 ! print *,'Before CALL RADIA0'
301 jscott 1.1 CALL RADIA0 (IO,JM,CO2,READGHG) 5090.
302 jscott 1.3 ! print *,'After CALL RADIA0'
303 jscott 1.1 INCHM=NRAD/NDYN 5091.
304     C**** CLOUD LAYER INDICES USED FOR DIAGNOSTICS 5092.
305     DO 43 L=1,LM 5093.
306     LLOW=L 5094.
307     IF (.5*(PLE(L+1)+PLE(L+2)).LT.786.) GO TO 44 5095.
308     43 CONTINUE 5096.
309     44 LMID1=LLOW+1 5097.
310     DO 45 L=LMID1,LM 5098.
311     LMID=L 5099.
312     IF (.5*(PLE(L+1)+PLE(L+2)).LT.430.) GO TO 46 5100.
313     45 CONTINUE 5101.
314     46 LHI1=LMID+1 5102.
315     LHI=LM 5103.
316     IF (LHI1.GT.LHI) LHI=LHI1 5104.
317     WRITE (6,47) LLOW,LMID1,LMID,LHI1,LHI 5105.
318     47 FORMAT (' LOW CLOUDS IN LAYERS 1-',I2,' MID LEVEL CLOUDS IN',5106.
319     * ' LAYERS',I3,'-',I2,' HIGH CLOUDS IN LAYERS',I3,'-',I2) 5107.
320     C**** NO RADIATION AVERAGING IJRA=1 JRA=1 IRA=1 5108.
321     C**** RADIATION AVERAGING IN I 2 1 2 5109.
322     C**** RADIATION AVERAGING IN I AND J 4 2 2 5110.
323     JRA=(IJRA+2)/3 5111.
324     IRA=IJRA/JRA 5112.
325     50 JALTER=MOD(NSTEP,NRAD*JRA)/NRAD 5113.
326     JDAYR=JDAY
327     JYEARR=JYEAR
328     IALTER=MOD(NSTEP,NRAD*IJRA)/(NRAD*JRA) 5114.
329     TNOW=JYEAR+(JDAY-.5)/365. 5127.1
330     IF (JDAY.NE.JDLAST.AND.KSTREND.NE.0) then
331     IF (KSTREND.GT.0) then
332     if(abs(CFS0X-1.0).gt.1.e3)then
333     print *,'Wrong setting KSTREND=',KSTREND,' CFS0X=',
334     & CFS0X
335     endif
336     S0X=S0X0*(1.+S0RATE/70.*(TNOW-INYEAR))
337     c print *,'From radia S0X0=',S0X0,' S0X=',S0X
338     c print *,'S0X=',S0X,' TNOW=',TNOW
339     ELSE
340 jscott 1.3 print *,'From radia TNOW=',TNOW
341 jscott 1.1 call obssolar(S0X,TNOW)
342 jscott 1.3 S0X=S0X/1367.
343 jscott 1.1 ENDIF
344     ENDIF
345     S0=S0X*1367./RSDIST 5115.
346     S0R=S0
347 jscott 1.3 ! print *,'S0=',S0
348 jscott 1.1 C**** CALCULATE AVERAGE COSINE OF ZENITH ANGLE FOR CURRENT COMP3 STEP 5116.
349     C**** AND RADIATION PERIOD 5117.
350     ROT1=TWOPI*TOFDAY/24. 5118.
351     if(DC25)then
352     ROT2=ROT1+TWOPI*DTCNDS/SDAY 5119.
353     CALL COSZT (IO,JM,SIND,COSD,ROT1,ROT2,COSZ1) 5120.
354     else
355     ROT2=ROT1+TWOPI
356     CALL COSZR (IO,JM,SIND,COSD,ROT1,ROT2,COSZ1)
357     end if
358     if(HPRNT)then
359     print *,' radia TAU=',TAU
360     print *,' CLDSS'
361     print *,(CLDSS(1,7,L),L=1,LM)
362     print *,' CLDMC'
363     print *,(CLDMC(1,7,L),L=1,LM)
364     endif
365     cprint *,' form radia TAU=',TAU,'MODRD=',MODRD
366     C
367     IF(MODRD.NE.0) GO TO 840 5121.
368     C
369     ROT2=ROT1+TWOPI*NRAD*DT/SDAY 5122.
370     CALL COSZS (IO,JM,SIND,COSD,ROT1,ROT2,COSZ2,COSZA) 5123.
371     C**** 5124.
372     C**** COMPUTE EARTH ALBEDOS AND OTHER PARAMETERS FOR BEGINNING OF DAY 5125.
373    
374     TNOW=JYEAR+(JDAY-.5)/365. 5127.1
375     KWRITE=0
376     if(JMONTH.ne.JMLAST) then
377     KWRITE=1
378     print *,' FROM radia READGHG=',READGHG
379     if(READGHG.eq.2) call tgases(CO2,JMONTH)
380     c print *,(zco2(1,j,1)
381     c & *28.97296245/44.0*1.e-3,j=1,jm)
382     !ppb(m) to kg per volume base
383     if(READGHG.eq.1) call rtgases(CO2,JMONTH)
384     #if ( defined PREDICTED_BC)
385 jscott 1.3 if(READGHG.eq.0)then
386 jscott 1.1 nrdbc=nrdbc+1
387     do l=1,lm
388     do j=1,jm
389     bcod(1,j,l)=bcodmn(j,l,nrdbc)
390     enddo
391     enddo
392     if(nrdbc.eq.12)then
393     nrdbc=0.
394     endif
395 jscott 1.3 endif
396 jscott 1.1 #endif
397     endif
398     JMLAST=JMONTH
399     IF (JDAY.NE.JDLAST.AND.KTREND.GT.0)
400     & CALL FORGET(TNOW,KTREND,KWRITE)
401     IF (JDAY.NE.JDLAST)then
402     #ifdef PREDICTED_GASES
403     call chemglobal(P)
404     #endif
405     #if ( defined PREDICTED_BC)
406     call sulfr(BSO4LAND,BSO4OCEAN,1986.)
407     do j=1,jm
408     FLAND=FDATA(1,J,2)
409     BSO4TOTAL(j)=BSO4LAND(j)*FLAND+BSO4OCEAN(j)*(1.-FLAND)
410     enddo
411     #endif
412     CALL RCOMPT
413     if(CLDFEED)then
414     DTSURFAV=0.
415     C
416     do j=1,jm
417 jscott 1.3 DTSURFAV=DTSURFAV+DT2MGL(J)*DXYP(j)
418 jscott 1.1 end do !j
419     DTSURFAV=DTSURFAV/AREAG
420     do j=1,jm
421     do k=1,3
422     cfcld(j,k)=1.+coefcl(k)*DTSURFAV
423     end do ! k
424     end do ! j
425     endif
426     ENDIF
427     JDLAST=JDAY 5129.
428     IHOUR=1.5+TOFDAY 5130.
429     if(CLDFEED)then
430     if (KWRITE.eq.1)then
431     print *,'cfcld'
432     print 9456,cfcld
433     print *,' DTSURF'
434 jscott 1.3 print 9456,DT2MGL
435 jscott 1.1 print *,' DTSURFAV=',DTSURFAV
436     9456 format(12f6.2)
437     endif
438     C Save clouds predicted by model
439     do k=1,LM
440     do j=1,JM
441     CLDSSF(j,k)=CLDSS(1,j,k)
442     CLDMCF(j,k)=CLDMC(1,j,k)
443     enddo
444     enddo
445     do k=1,LM
446     if(k.le.5)then
447     k1=1
448     else
449     k1=2
450     endif
451     do j=1,JM
452     CLDSS(1,j,k)=cfcld(j,k1)*CLDSS(1,j,k)
453     CLDMC(1,j,k)=cfcld(j,3)*CLDMC(1,j,k)
454     enddo
455     enddo
456     endif
457     C
458     CB READING OF CLOUD
459     if(abs(PCLOUD-3.).lt.1.5)then
460     910 continue
461     if(nreadcld.eq.nrcldmax)go to 900
462     read(585,END=900)TFDAYF,JDATEF,JMNTHF,CLDSSF,CLDMCF,IRAND
463     nreadcld=nreadcld+1
464     if(IFIRST.eq.1)then
465     print *,' radia.f PCLOUD=',PCLOUD
466     if(PCLOUD.eq.2)print *,' FIXED MC and SS CLOUDS'
467     if(PCLOUD.eq.4)print *,' FIXED MC CLOUDS ONLY'
468     if(PCLOUD.eq.3)print *,' FIXED SS CLOUDS ONLY'
469     print *,TOFDAY,JDATE,JMONTH
470     print *,TFDAYF,JDATEF,JMNTHF
471     print *,' DTCNDS=',DTCNDS/3600.
472     print *,' DT*NRAD=',DT*NRAD/3600.
473     if(.not.WRCLD)IFIRST=0
474     endif
475     if(abs(TOFDAY-TFDAYF).gt.1.e-3.or.JDATE.ne.JDATEF.or.
476     * JMONTH.ne.JMNTHF)then
477     print *,' RADIA, disagrement in clouds'
478     print *,TOFDAY,JDATE,JMONTH
479     print *,TFDAYF,JDATEF,JMNTHF
480     stop
481     endif
482     go to 920
483     900 rewind 585
484     nreadcld=0
485     print *,' END OF file85'
486     print *,JYEAR
487     print *,TOFDAY,JDATE,JMONTH
488     print *,' REWIND 85'
489     go to 910
490     920 continue
491     c if (KWRITE.eq.1)then
492     c print *,'cfcld'
493     c print 9456,cfcld
494     c print *,' DTSURF'
495     c print 9456,DTSURF
496     c9456 format(12f6.2)
497     c endif
498     CALL RINIT (IRAND)
499     do 930 k=1,LM
500     do 930 j=1,JM
501     if(PCLOUD.ne.4)CLDSS(1,j,k)=CLDSSF(j,k)
502     if(PCLOUD.ne.3)CLDMC(1,j,k)=CLDMCF(j,k)
503     930 continue
504     endif
505     CE END OF READING OF CLOUD
506     if(WRCLD)then
507     if(NWRCLD.eq.1)then
508     CALL RFINAL(IRAND)
509     if(IFIRST.eq.1)print *,' SHORT CLOUDS RECORD'
510     write(581)TOFDAY,JDATE,JMONTH,CLDSS,CLDMC,IRAND
511     elseif(NWRCLD.eq.2)then
512     if(IFIRST.eq.1)print *,' LONG CLOUDS RECORD'
513     do 1150 k=1,14
514     do 1150 j=1,JM0
515     if(k.le.2)then
516     ODATA2(j,k)=ODATA(1,j,k)
517     BDATA2(j,k)=BLDATA(1,j,k)
518     FDATA2(j,k)=FDATA(1,j,k+1)
519     endif
520     if(k.le.3)RQT2(j,k)=RQT(1,j,k)
521     GDATA2(j,k)=GDATA(1,j,k)
522     1150 continue
523     CALL RFINAL(IRAND)
524     write(581)TOFDAY,JDATE,JMONTH,CLDSS,CLDMC,IRAND,
525     * JDAY,JYEAR,T,Q,P,
526     * ODATA2,BDATA2,FDATA2,GDATA2,RQT2
527     & ,CLDSSF,CLDMCF
528 jscott 1.3 #if ( defined CLM )
529     write(581)asdirclm,asdifclm,aldirclm,aldifclm
530     #endif
531    
532 jscott 1.1 else
533     print *,' NWRCLD=',NWRCLD
534     stop
535     endif
536     IFIRST=0
537     endif
538     C**** 5131.
539     C**** MAIN J LOOP 5132.
540     C**** 5133.
541     DO 600 J=1,JM 5134.
542     IF ((J-1)*(JM-J).NE.0) GO TO 140 5135.
543     C**** CONDITIONS AT THE POLES 5136.
544     POLE=.TRUE. 5137.
545     MODRJ=0 5138.
546     IMAX=1 5139.
547     GO TO 160 5140.
548     C**** CONDITIONS AT NON-POLAR POINTS 5141.
549     140 POLE=.FALSE. 5142.
550     MODRJ=MOD(J+JALTER,JRA) 5143.
551     IMAX=IM 5144.
552     160 XFRADJ=.2+1.2*COSP(J)*COSP(J) 5145.
553     JLAT=J 5145.1
554     if(GSOEQ)then
555     C
556     RVOL=0.012
557     FVOL=0.0045
558     FGOLDU(1)=(RVOL+FVOL)/RVOL
559     C
560     endif
561     IF(MODRJ.EQ.0) CALL RCOMPJ 5146.
562 jscott 1.3 SWIN(j,1)=0.0
563     SWNET(j,1)=0.0
564     SWIN(j,2)=0.0
565     SWNET(j,2)=0.0
566 jscott 1.1 C**** 5147.
567     C**** MAIN I LOOP 5148.
568     C**** 5149.
569     IM1=IM 5150.
570     DO 500 I=1,IMAX 5151.
571     MODRIJ=MODRJ+MOD(I+IALTER,IRA) 5152.
572     IF(POLE) MODRIJ=0 5153.
573     JR=J
574     C**** DETERMINE FRACTIONS FOR SURFACE TYPES AND COLUMN PRESSURE 5155.
575     PLAND=FDATA(I,J,2) 5156.
576     PWATER=1.-PLAND
577     POICE=ODATA(I,J,2)*(1.-PLAND) 5157.
578     POCEAN=(1.-PLAND)-POICE 5158.
579     if(POCEAN.LE.1.E-5)then
580     POCEAN=0.
581     POICE=PWATER
582     endif
583     PLICE=FDATA(I,J,3)*PLAND 5159.
584     PEARTH=PLAND-PLICE 5160.
585     SP=P(I,J) 5161.
586     C**** 5162.
587     C**** DETERMINE CLOUDS (AND THEIR OPTICAL DEPTHS) SEEN BY RADIATION 5163.
588     C**** 5164.
589     X=999999. 5164.1
590     c RANDSS=RANDU(X) 5165.
591     c RANDMC=RANDU(X) 5166.
592     CALL RANDUU(RANDSS,X)
593     CALL RANDUU(RANDMC,X)
594     if(HPRNT)then
595     print *,' radia TAU=',TAU
596     print *,' RANDSS=',RANDSS
597     print *,' RANDMC=',RANDMC
598     endif
599     C
600     CSS=0. 5167.
601     CMC=0. 5168.
602     DEPTH=0. 5169.
603     LTOP=0 5169.1
604     DO 210 L=1,LM 5170.
605     RTAU(L)=0. 5171.
606     210 TOTCLD(L)=0. 5172.
607     DO 240 L=1,LM 5173.
608     IF(CLDSS(I,J,L).LT.RANDSS) GO TO 220 5174.
609     RTAUSS=.013333*(PTOP-100.+SIG(L)*SP) 5175.
610     IF(RTAUSS.LT.0.) RTAUSS=0. 5176.
611     IF (T(I,J,L)*PK(I,J,L).LT.TCIR) RTAUSS=.3333333 5177.
612     RTAU(L)=RTAUSS 5178.
613     CSS=1. 5179.
614     AJL(J,L,28)=AJL(J,L,28)+CSS 5180.
615     TOTCLD(L)=1. 5181.
616     LTOP=L 5181.1
617     220 IF(CLDMC(I,J,L).LE.RANDMC) GO TO 240 5182.
618     RTAUMC=DSIG(L)*SP*.08 5183.
619     IF(RTAUMC.GT.RTAU(L)) RTAU(L)=RTAUMC 5184.
620     CMC=1. 5185.
621     AJL(J,L,29)=AJL(J,L,29)+CMC 5186.
622     TOTCLD(L)=1. 5187.
623     LTOP=L 5187.1
624     DEPTH=DEPTH+SP*DSIG(L) 5188.
625     AJL(J,L,19)=AJL(J,L,19)+TOTCLD(L) 5189.
626     240 CONTINUE
627     AJ(J,57)=AJ(J,57)+CSS*POCEAN 5190.
628     BJ(J,57)=BJ(J,57)+CSS*PLAND 5191.
629     CJ(J,57)=CJ(J,57)+CSS*POICE 5192.
630     DJ(JR,57)=DJ(JR,57)+CSS*DXYP(J) 5193.
631     AJ(J,58)=AJ(J,58)+CMC*POCEAN 5194.
632     BJ(J,58)=BJ(J,58)+CMC*PLAND 5195.
633     CJ(J,58)=CJ(J,58)+CMC*POICE 5196.
634     DJ(JR,58)=DJ(JR,58)+CMC*DXYP(J) 5197.
635     AIJ(I,J,17)=AIJ(I,J,17)+CMC 5198.
636     AJ(J,80)=AJ(J,80)+DEPTH*POCEAN 5199.
637     BJ(J,80)=BJ(J,80)+DEPTH*PLAND 5200.
638     CJ(J,80)=CJ(J,80)+DEPTH*POICE 5201.
639     DJ(JR,80)=DJ(JR,80)+DEPTH*DXYP(J) 5202.
640     CLDCV=CMC+CSS-CMC*CSS 5203.
641     cldd4tem(j)=cldd4tem(j)+CLDCV
642     ncldd4tem(j)=ncldd4tem(j)+1
643     AJ(J,59)=AJ(J,59)+CLDCV*POCEAN 5204.
644     BJ(J,59)=BJ(J,59)+CLDCV*PLAND 5205.
645     CJ(J,59)=CJ(J,59)+CLDCV*POICE 5206.
646     DJ(JR,59)=DJ(JR,59)+CLDCV*DXYP(J) 5207.
647     C**** clear sky condinion
648     CSS0=CSS
649     CMC0=CMC
650     do L=1,LM
651     RTAU0(L)=RTAU(L)
652     TOTCLD0(L)=TOTCLD(L)
653     enddo
654     if(CMC.le.0.and.CSS.le.0)then
655     CLEAR(J)=1
656     else
657     CLEAR(J)=0
658     endif
659     if(STRARFOR.or.CO2FOR.or.S0FOR.or.FORBC)then
660     CLEAR(J)=0
661     endif
662     C**** 5247.
663     C**** SET UP VERTICAL ARRAYS OMITTING THE I AND J INDICES 5248.
664     C**** 5249.
665     C**** EVEN PRESSURES 5250.
666     DO 340 L=1,LM 5251.
667     PLE(L)=SIGE(L)*SP+PTOP 5252.
668     C**** TEMPERATURES 5253.
669     TL(L)=T(I,J,L)*PK(I,J,L) 5254.
670     C**** MOISTURE VARIABLES 5255.
671     QL(L)=Q(I,J,L) 5256.
672     340 CONTINUE 5257.
673     C**** 5258.
674     C**** RADIATION, SOLAR AND THERMAL 5259.
675     C**** 5260.
676    
677     DO 420 K=1,3 5261.
678     420 TL(LM+K)=RQT(I,J,K) 5262.
679     COSZ=COSZA(I,J) 5263.
680     TGO=ODATA(I,J,1)+TF 5264.
681     TGOI=GDATA(I,J,3)+TF 5265.
682     TGLI=GDATA(I,J,13)+TF 5266.
683     TGE=GDATA(I,J,4)+TF 5267.
684     TS=BLDATA(I,J,2) 5268.
685     SNOWOI=GDATA(I,J,1) 5269.
686     SNOWLI=GDATA(I,J,12) 5270.
687     SNOWE=GDATA(I,J,2) 5271.
688     AGESN=GDATA(I,J,11) 5272.
689     c print *,'From radia J=',J,' TAU=',TAU,' TS=',TS
690     c print *,TGO,TGOI,TGLI,TGE
691     c print *,SNOWOI,SNOWLI,SNOWE,AGESN
692     WEARTH=(GDATA(I,J,5)+GDATA(I,J,6))/(VDATA(I,J,9)+1.E-20) 5273.
693     DO 430 K=1,8 5274.
694     430 PVT(K)=VDATA(I,J,K) 5275.
695     WS=BLDATA(I,J,1) 5276.
696     do 439 L=1,LM+1
697     SRHR(I,J,L)=0.
698     TRHR(I,J,L)=0.
699     if(L.le.4)then
700     SNFS(I,J,L)=0.
701     TNFS(I,J,L)=0.
702     if(L.le.3)then
703     SRHRS(I,J,L)=0.
704     TRHRS(I,J,L)=0.
705     endif
706     endif
707     439 continue
708     TRNFP0(J)=0.
709     TRNFP1(J)=0.
710     TRINCG(I,J)=0.
711     BTMPW(I,J)=0.
712     SRDAN=0.
713     SRNAN=0.
714     do 449 K=1,9
715     ALB(I,J,K)=0.
716     ALBJ(J,K)=0.
717     449 continue
718     do 499 ii=1,3
719     COSZ=COSZA(I,J)
720     PLAND=FDATA(I,J,2)
721     PWATER=1.-PLAND
722     POICE=ODATA(I,J,2)*(1.-PLAND)
723     POCEAN=(1.-PLAND)-POICE
724     if(POCEAN.LE.1.E-5)then
725     POCEAN=0.
726     POICE=PWATER
727     endif
728     PLICE=FDATA(I,J,3)*PLAND
729     PEARTH=PLAND-PLICE
730     if(ii.eq.1)then
731     BSO4BC=BSO4OCEAN(J)/BSO4TOTAL(J)
732     PTYPE=POCEAN
733     POICE=0.
734     POCEAN=1.
735     PLAND=0.
736     PEARTH=0.
737     PLICE=0.
738     TGAL=0.
739     else if(ii.eq.3)then
740     BSO4BC=BSO4OCEAN(J)/BSO4TOTAL(J)
741     PTYPE=POICE
742     POICE=1.
743     POCEAN=0.
744     PLAND=0.
745     PEARTH=0.
746     PLICE=0.
747     TGAL=TGOI
748     else
749     BSO4BC=BSO4LAND(J)/BSO4TOTAL(J)
750     PTYPE=PLAND
751     POCEAN=0.
752     POICE=0.
753     PWATER=0.
754     PLICE=FDATA(I,J,3)
755     PEARTH=1.-PLICE
756     TGAL=TGE*PEARTH+TGLI*PLICE
757     PLAND=1.
758     endif
759     if(PTYPE.lt.1.e-10)go to 499
760     if(ii.gt.1)then
761     if(TGAL.lt.263.)then
762     FRSNALB=0.30
763     elseif(TGAL.lt.273.)then
764     FRSNALB=0.30-0.015*(TGAL-263.)
765     else
766     FRSNALB=0.15
767     endif
768     endif !ii
769     c FGOLDU(2)=XFRADJ*(1.-PEARTH) 5277.
770     c FGOLDU(3)=XFRADJ*PEARTH 5278.
771     FGOLDU(2)=XFRADJ*(1.-PLAND)
772     FGOLDU(3)=XFRADJ*PLAND
773     if(CO2FOR)then
774     IF (KTREND.GT.0) THEN
775     CALL FORGET(TNOW,KTREND,0)
776     ELSE
777     FULGAS(2)=CO2
778     ENDIF
779     endif
780     !
781 jscott 1.3 #if ( defined IPCC_EMI )
782     FULGAS(2)=CO2
783     #endif
784    
785 jscott 1.1 #if ( defined FIXED_FOR )
786     FULGAS(2)=1.
787     #endif
788     !
789     if(S0FOR)then
790     S0=S0R
791     endif
792     CSS=CSS0
793     CMC=CMC0
794     do L=1,LM
795     RTAU(L)=RTAU0(L)
796     TOTCLD(L)=TOTCLD0(L)
797     enddo
798    
799     Chemstry Model Patch 092295
800     c
801     PTYPER=max(PTYPE,5.e-3)
802     FAERSOL=PLAND/PTYPER*3.0*CFAEROSOL
803     C Factor 3 is a "parameterization" for indirec aerosol effect
804     c For normal setting CFAEROSOL=1.0 (default)
805     c For hight forcing CFAEROSOL=0.5
806     c For low forcing CFAEROSOL=2.0
807     c
808     c ==========================
809     #if ( defined PREDICTED_BC)
810     FBC=BSO4BC*CFBC
811     #endif
812     ILON=I 5278.1
813     JLAT=J 5278.2
814     if(GHSFALB)then
815     BVSURFA=ALBCF
816     XVSURFA=ALBCF
817     BNSURFA=ALBCF
818     XNSURFA=ALBCF
819     else
820     BVSURFA=0.0
821     XVSURFA=0.0
822     BNSURFA=0.0
823     XNSURFA=0.0
824     endif
825     if(STRARFOR)then
826     c if(j.gt.JM/2)then
827     FGOLDU(1)=FGOLDU1
828     c else
829     c FGOLDU(1)=1.
830     c endif
831     endif
832     if(J.eq.-10)then
833     print *,'From radia GHSFALB=',GHSFALB,' ALBCF=',ALBCF
834     print *,'From radia BVSURFA=',BVSURFA
835     endif
836     c if(J.eq.1)then
837     c print *,ii,PTYPE,PWATER,POCEAN,POICE,PLAND,PEARTH,PLICE
838     c print *,(FGOLDU(if),if=1,5)
839     c endif
840     c print *,'S0 for real calculation=',S0
841     CALL RCOMPX 5279.
842     if (IRFIRST.eq.1.and.READGHG.eq.1)then
843     CALL WRITER(12)
844     if(ii.ge.2)IRFIRST=0
845     endif
846     IF(DMOD(TAU,365.*24.).EQ.0..and.J.eq.JM/2) then
847     print *,' tau=',TAU,' J=',J
848     CALL WRITER (1,0)
849     endif
850     c CALL WRITER (13) 5280.
851     SRHR(I,J,1)=SRHR(I,J,1)+SRNFLB(1)*PTYPE
852     c GHOST FORCINGS
853     if(GHSF)then
854    
855     do LFR=1,lm+1
856     if(LFR.gt.1)then
857     TRFCRL(LFR-1)=TRFCRL(LFR-1)-ghostf(LFR,J)
858     TRNFLB(LFR)=TRNFLB(LFR)+ghflux(LFR,J)
859     endif
860     enddo
861     endif
862     c
863     TRHR(I,J,1)=TRHR(I,J,1)+(STBO*(POCEAN*TGO**4+POICE*TGOI**4
864     * +PLICE*TGLI**4+PEARTH*TGE**4)-TRNFLB(1))*PTYPE
865     if(GHSF)then
866     TRHR(I,J,1)=TRHR(I,J,1)+ghostf(1,J)*PTYPE
867     TRDFLB(1)=TRDFLB(1)+ghostf(1,J)
868     endif
869     C *****
870     TRSURF(J,ii)=STBO*(POCEAN*TGO**4+POICE*TGOI**4
871     * +PLICE*TGLI**4+PEARTH*TGE**4)-TRNFLB(1)
872 jscott 1.3 ! if(TAU.gt.3623.0.and.TAU.lt.3744.0) then
873     if(ii.eq.-2) then
874     if(J.eq.29) then
875     ! print *,'From radia TAU=',TAU
876     ! print *,'J=',j,' ii=',ii
877     ! print *,POCEAN,POICE,PLICE,PEARTH
878     ! print *,TGO,TGOI,TGLI,TGE
879     ! print *,' TRNFLB(1)=',TRNFLB(1)
880     ! print *,'TRSURF(J,ii)=',TRSURF(J,ii),' TRDFLB(1)=',TRDFLB(1)
881     ! print *,' TRUFLB(1)=',TRUFLB(1),' SRBO*T**4=',
882     ! & STBO*(POCEAN*TGO**4+POICE*TGOI**4
883     ! * +PLICE*TGLI**4+PEARTH*TGE**4)
884     ! print *,'CLD',CLDCV,CMC,CSS
885     ! if(ii.eq.2)print *,'TGL=',TGE,STBO*TGE**4
886     ! write(99,*),'From radia ',TAU,j,ii
887     ! write(99,*), CLDCV,CMC,CSS
888     write(129), TAU,TOFDAY,CLDCV,TGE,TRSURF(J,ii),SRDFLB(1)
889     ! write(99,*), TAU,TRSURF(J,ii),
890     ! & STBO*(POCEAN*TGO**4+POICE*TGOI**4
891     ! * +PLICE*TGLI**4+PEARTH*TGE**4),
892     ! & TRNFLB(1)
893     endif
894     if(J.eq.30) then
895     write(130), TAU,TOFDAY,CLDCV,TGE,TRSURF(J,ii),SRDFLB(1)
896     endif
897     endif
898     ! endif
899 jscott 1.1 if(GHSF)then
900     TRSURF(J,ii)=TRSURF(J,ii)+ghostf(1,J)
901     endif
902     SRSURF(J,ii)=SRNFLB(1)
903     DO 440 L=1,LM 5284.
904     SRHR(I,J,L+1)=SRHR(I,J,L+1)+SRFHRL(L)*PTYPE
905     440 TRHR(I,J,L+1)=TRHR(I,J,L+1)-TRFCRL(L)*PTYPE
906     DO 450 LR=1,3 5287.
907     SRHRS(I,J,LR)=SRHRS(I,J,LR)+SRFHRL(LM+LR)*PTYPE
908     450 TRHRS(I,J,LR)=TRHRS(I,J,LR)-TRFCRL(LM+LR)*PTYPE
909     DO 460 K=1,4 5290.
910     SNFS(I,J,K)=SNFS(I,J,K)+SRNFLB(K+LM)*PTYPE
911     460 TNFS(I,J,K)=TNFS(I,J,K)+(TRNFLB(K+LM)-TRNFLB(1))*PTYPE
912     TRNFP0(J)=TRNFP0(J)+TRNFLB(4+LM)*PTYPE
913     TRNFP1(J)=TRNFP1(J)+TRNFLB(1+LM)*PTYPE
914     TRINCG(I,J)=TRINCG(I,J)+TRDFLB(1)*PTYPE
915     BTMPW(I,J)=BTMPW(I,J)+(BTEMPW-TF)*PTYPE
916     c ALB(I,J,1)=SRNFLB(1)/(SRDFLB(1)+1.E-20) 5295.
917     SRDAN=SRDAN+SRDFLB(1)*PTYPE
918     SRNAN=SRNAN+SRNFLB(1)*PTYPE
919     ALB(I,J,2)=ALB(I,J,2)+PLAVIS*PTYPE
920     ALB(I,J,3)=ALB(I,J,3)+PLANIR*PTYPE
921     ALB(I,J,4)=ALB(I,J,4)+ALBVIS*PTYPE
922     ALB(I,J,5)=ALB(I,J,5)+ALBNIR*PTYPE
923     ALB(I,J,6)=ALB(I,J,6)+SRRVIS*PTYPE
924     ALB(I,J,7)=ALB(I,J,7)+SRRNIR*PTYPE
925     ALB(I,J,8)=ALB(I,J,8)+SRAVIS*PTYPE
926     ALB(I,J,9)=ALB(I,J,9)+SRANIR*PTYPE
927     ALB1=SRNFLB(1)/(SRDFLB(1)+1.E-20)
928     C **********
929     ALBJ(J,2)=PLAVIS
930     ALBJ(J,3)=PLANIR
931     ALBJ(J,4)=ALBVIS
932     ALBJ(J,5)=ALBNIR
933     ALBJ(J,6)=SRRVIS
934     ALBJ(J,7)=SRRNIR
935     ALBJ(J,8)=SRAVIS
936     ALBJ(J,9)=SRANIR
937     ALBJ(J,1)=SRNFLB(1)/(SRDFLB(1)+1.E-20)
938     C *********
939     COSZ=COSZ2(I,J)
940     do L=LM,1,-1
941     AJL(J,L,43)=AJL(J,L,43)-TRNFLB(L+1)*PTYPE
942     AJL(J,L,42)=AJL(J,L,42)+(SRNFLB(L+1)*COSZ)*PTYPE
943     enddo
944     if(ii.eq.2)then
945     #if ( defined CLM )
946     C for TEM CLM
947     DSWSRF(j)=SRDFLB(1)
948     DLWSRF(j)=TRSURF(J,2)
949 jscott 1.3 ULWSRF(j)=TRUFLB(1)
950 jscott 1.1 DSWVIS(j)=SRDVIS
951     DSWNIR(j)=SRDNIR
952     C for TEM CLM
953     #endif
954     PLAND=PTYPE
955     BJ(J,1)=BJ(J,1)+(S0*COSZ)*PLAND
956     BJ(J,2)=BJ(J,2)+(SRNFLB(4+LM)*COSZ)*PLAND
957     c BSNFS1=BSNFS1+(SRNFLB(1+LM)*COSZ)*PLAND
958     c BJ(J,5)=BJ(J,5)+(SRNFLB(1)*COSZ/(ALB1+1.E-20))*PLAND
959     BJ(J,5)=BJ(J,5)+(SRDFLB(1)*COSZ)*PLAND
960     BJ(J,6)=BJ(J,6)+(SRNFLB(1)*COSZ)*PLAND
961     BJ(J,55)=BJ(J,55)+(BTEMPW-TF)*PLAND
962     BJ(J,67)=BJ(J,67)+TRDFLB(1)*PLAND
963     BJ(J,70)=BJ(J,70)-(TRNFLB(4+LM)-TRNFLB(1))*PLAND
964     BJ(J,7)=BJ(J,7)-TRNFLB(4+LM)*PLAND
965     BJ(J,8)=BJ(J,8)-TRNFLB(1+LM)*PLAND
966     BJ(J,3)=BJ(J,3)+(SRNFLB(1+LM)*COSZ)*PLAND
967     BJ(J,71)=BJ(J,71)-(TRNFLB(1+LM)-TRNFLB(1))*PLAND
968     DO 761 K=2,9
969     BJ(J,K+70)=BJ(J,K+70)+(S0*COSZ)*ALBJ(J,K)*PLAND
970     761 CONTINUE
971     else if(ii.eq.1)then
972     POCEAN=PTYPE
973     AJ(J,1)=AJ(J,1)+(S0*COSZ)*POCEAN
974     AJ(J,2)=AJ(J,2)+(SRNFLB(4+LM)*COSZ)*POCEAN
975     AJ(J,5)=AJ(J,5)+(SRDFLB(1)*COSZ)*POCEAN
976     AJ(J,6)=AJ(J,6)+(SRNFLB(1)*COSZ)*POCEAN
977     AJ(J,55)=AJ(J,55)+(BTEMPW-TF)*POCEAN
978     AJ(J,67)=AJ(J,67)+TRDFLB(1)*POCEAN
979     AJ(J,70)=AJ(J,70)-(TRNFLB(4+LM)-TRNFLB(1))*POCEAN
980     AJ(J,7)=AJ(J,7)-TRNFLB(4+LM)*POCEAN
981     AJ(J,8)=AJ(J,8)-TRNFLB(1+LM)*POCEAN
982     AJ(J,3)=AJ(J,3)+(SRNFLB(1+LM)*COSZ)*POCEAN
983     AJ(J,71)=AJ(J,71)-(TRNFLB(1+LM)-TRNFLB(1))*POCEAN
984     #if ( defined OCEAN_3D )
985 jscott 1.3 SWIN(j,1)=SRDFLB(1)
986     SWNET(j,1)=SRNFLB(1)
987 jscott 1.1 #endif
988     C
989     DO K=2,9
990     AJ(J,K+70)=AJ(J,K+70)+(S0*COSZ)*ALBJ(J,K)*POCEAN
991     END DO
992     C
993     else
994     POICE=PTYPE
995     CJ(J,1)=CJ(J,1)+(S0*COSZ)*POICE
996     CJ(J,2)=CJ(J,2)+(SRNFLB(4+LM)*COSZ)*POICE
997     CJ(J,5)=CJ(J,5)+(SRDFLB(1)*COSZ)*POICE
998     CJ(J,6)=CJ(J,6)+(SRNFLB(1)*COSZ)*POICE
999     CJ(J,55)=CJ(J,55)+(BTEMPW-TF)*POICE
1000     CJ(J,67)=CJ(J,67)+TRDFLB(1)*POICE
1001     CJ(J,70)=CJ(J,70)-(TRNFLB(4+LM)-TRNFLB(1))*POICE
1002     CJ(J,7)=CJ(J,7)-TRNFLB(4+LM)*POICE
1003     CJ(J,8)=CJ(J,8)-TRNFLB(1+LM)*POICE
1004     CJ(J,3)=CJ(J,3)+(SRNFLB(1+LM)*COSZ)*POICE
1005     CJ(J,71)=CJ(J,71)-(TRNFLB(1+LM)-TRNFLB(1))*POICE
1006     #if ( defined OCEAN_3D )
1007 jscott 1.3 SWIN(j,2)=SRDFLB(1)
1008     SWNET(j,2)=SRNFLB(1)
1009 jscott 1.1 #endif
1010     C
1011     DO K=2,9
1012     CJ(J,K+70)=CJ(J,K+70)+(S0*COSZ)*ALBJ(J,K)*POICE
1013     END DO
1014     endif
1015     if(CLEAR(J).eq.0)then
1016 jscott 1.3 PLAND=FDATA(I,J,2)
1017     PWATER=1.-PLAND
1018     POICE=ODATA(I,J,2)*(1.-PLAND)
1019     POCEAN=(1.-PLAND)-POICE
1020     if(POCEAN.LE.1.E-5)then
1021     POCEAN=0.
1022     POICE=PWATER
1023     endif
1024     PLICE=FDATA(I,J,3)*PLAND
1025     PEARTH=PLAND-PLICE
1026     if(ii.eq.1)then
1027     PTYPE=POCEAN
1028     POICE=0.
1029     POCEAN=1.
1030     PLAND=0.
1031     PEARTH=0.
1032     PLICE=0.
1033     else if(ii.eq.3)then
1034     PTYPE=POICE
1035     POICE=1.
1036     POCEAN=0.
1037     PLAND=0.
1038     PEARTH=0.
1039     PLICE=0.
1040     else
1041     PTYPE=PLAND
1042     POCEAN=0.
1043     POICE=0.
1044     PWATER=0.
1045     PLICE=FDATA(I,J,3)
1046     PEARTH=1.-PLICE
1047     PLAND=1.
1048     endif
1049     COSZ=COSZA(I,J)
1050 jscott 1.1 if(STRARFOR)then
1051     FGOLDU(1)=1.0
1052     elseif(CO2FOR)then
1053     FULGAS(2)=CO2F
1054     elseif(S0FOR)then
1055     c print *,'For S0FOR ',S0X0,CFS0X,S0X
1056     S0=S0X0/CFS0X*1367./RSDIST
1057     c print *,'S0 for forcing calculation=',S0
1058     elseif(FORBC)then
1059     FBC=0.0
1060     else
1061     CSS=0.
1062     CMC=0.
1063     DEPTH=0.
1064     LTOP=0.
1065     do 1210 L=1,LM
1066     RTAU(L)=0.
1067     TOTCLD(L)=0.
1068     1210 continue
1069     endif
1070     568 continue
1071     CALL RCOMPX
1072 jscott 1.3 endif ! CLEAR
1073 jscott 1.1 SRHRCL(J)=SRNFLB(1)
1074     TRHRCL(J)=-TRNFLB(1)
1075     ALBCL(J)=SRNFLB(1)/(SRDFLB(1)+1.e-20)
1076     SNP1CL(J)=SRNFLB(LM+1)
1077     SNP0CL(J)=SRNFLB(LM+4)
1078     TRINCL(J)=TRDFLB(1)
1079     TRP0CL(J)=TRNFLB(LM+4)
1080     TRP1CL(J)=TRNFLB(LM+1)
1081 jscott 1.3 ! if(J.eq.29) then
1082     ! if(TAU.gt.3623.0.and.TAU.lt.3744.0) then
1083     ! print *,' TRHRCL(j)=',TRHRCL(j)
1084     ! print *,' TRINCL(j)=',TRINCL(j)
1085     ! endif
1086     ! endif
1087 jscott 1.1 COSZ=COSZ2(I,J)
1088     do L=LM,1,-1
1089     AJL(J,L,46)=AJL(J,L,46)-TRNFLB(L+1)*PTYPE
1090     AJL(J,L,45)=AJL(J,L,45)+(SRNFLB(L+1)*COSZ)*PTYPE
1091     enddo
1092     C *******
1093     NCLR(J)=NCLR(J)+1
1094     SRHE5=SRHRCL(J)*COSZ/(ALBCL(J)+1.E-20)
1095     if(ii.eq.1)then
1096     POCEAN=PTYPE
1097     AJCLR(J,1)=AJCLR(J,1)+(S0*COSZ)*POCEAN
1098     AJCLR(J,2)=AJCLR(J,2)+(SNP0CL(J)*COSZ)*POCEAN
1099     AJCLR(J,4)=AJCLR(J,4)+(SRHRCL(J)*COSZ)*POCEAN
1100     AJCLR(J,5)=AJCLR(J,5)+SRHE5*POCEAN
1101     AJCLR(J,6)=AJCLR(J,6)+TRINCL(J)*POCEAN
1102     AJCLR(J,8)=AJCLR(J,8)-TRP0CL(J)*POCEAN
1103     AJCLR(J,9)=AJCLR(J,9)-TRP1CL(J)*POCEAN
1104     AJCLR(J,3)=AJCLR(J,3)+(SNP1CL(J)*COSZ)*POCEAN
1105     AJCLR(J,7)=AJCLR(J,7)+TRHRCL(J)*POCEAN
1106     else if(ii.eq.2)then
1107     PLAND=PTYPE
1108     BJCLR(J,1)=BJCLR(J,1)+(S0*COSZ)*PLAND
1109     BJCLR(J,2)=BJCLR(J,2)+(SNP0CL(J)*COSZ)*PLAND
1110     BJCLR(J,4)=BJCLR(J,4)+(SRHRCL(J)*COSZ)*PLAND
1111     BJCLR(J,5)=BJCLR(J,5)+SRHE5*PLAND
1112     BJCLR(J,6)=BJCLR(J,6)+TRINCL(J)*PLAND
1113     BJCLR(J,8)=BJCLR(J,8)-TRP0CL(J)*PLAND
1114     BJCLR(J,9)=BJCLR(J,9)-TRP1CL(J)*PLAND
1115     BJCLR(J,3)=BJCLR(J,3)+(SNP1CL(J)*COSZ)*PLAND
1116     BJCLR(J,7)=BJCLR(J,7)+TRHRCL(J)*PLAND
1117     else
1118     POICE=PTYPE
1119     CJCLR(J,1)=CJCLR(J,1)+(S0*COSZ)*POICE
1120     CJCLR(J,2)=CJCLR(J,2)+(SNP0CL(J)*COSZ)*POICE
1121     CJCLR(J,4)=CJCLR(J,4)+(SRHRCL(J)*COSZ)*POICE
1122     CJCLR(J,5)=CJCLR(J,5)+SRHE5*POICE
1123     CJCLR(J,6)=CJCLR(J,6)+TRINCL(J)*POICE
1124     CJCLR(J,8)=CJCLR(J,8)-TRP0CL(J)*POICE
1125     CJCLR(J,9)=CJCLR(J,9)-TRP1CL(J)*POICE
1126     CJCLR(J,3)=CJCLR(J,3)+(SNP1CL(J)*COSZ)*POICE
1127     CJCLR(J,7)=CJCLR(J,7)+TRHRCL(J)*POICE
1128     endif
1129     C *********
1130     499 continue
1131     ALB(I,J,1)=SRNAN/(SRDAN+1.E-20)
1132     500 IM1=I 5304.
1133     C**** 5305.
1134     C**** END OF MAIN LOOP FOR I INDEX 5306.
1135     C**** 5307.
1136     600 CONTINUE 5345.
1137     C**** 5346.
1138     C**** END OF MAIN LOOP FOR J INDEX 5347.
1139     C**** 5348.
1140     C**** ACCUMULATE THE RADIATION DIAGNOSTICS 5394.
1141     C**** 5395.
1142     700 DO 780 J=1,JM 5396.
1143     DXYPJ=DXYP(J) 5397.
1144     IMAX=IM 5398.
1145     IF(J.EQ.1.OR.J.EQ.JM) IMAX=1 5399.
1146     DO 720 L=1,LM 5400.
1147     ASRHR=0. 5401.
1148     ATRHR=0. 5402.
1149     DO 710 I=1,IMAX 5403.
1150     ASRHR=ASRHR+SRHR(I,J,L+1)*COSZ2(I,J) 5404.
1151     710 ATRHR=ATRHR+TRHR(I,J,L+1) 5405.
1152     AJL(J,L,9)=AJL(J,L,9)+ASRHR 5406.
1153     720 AJL(J,L,10)=AJL(J,L,10)+ATRHR 5407.
1154     ASNFS1=0. 5408.
1155     BSNFS1=0. 5409.
1156     CSNFS1=0. 5410.
1157     ATNFS1=0. 5411.
1158     BTNFS1=0. 5412.
1159     CTNFS1=0. 5413.
1160     DO 770 I=1,IMAX 5414.
1161     SP=P(I,J) 5415.
1162     COSZ=COSZ2(I,J) 5416.
1163     PLAND=FDATA(I,J,2) 5417.
1164     PWATER=1.-PLAND
1165     POICE=ODATA(I,J,2)*(1.-PLAND) 5418.
1166     POCEAN=(1.-PLAND)-POICE 5419.
1167     if(POCEAN.LE.1.E-5)then
1168     POCEAN=0.
1169     POICE=PWATER
1170     endif
1171     JR=J
1172     DO 740 LR=1,3 5421.
1173     ASJL(J,LR,3)=ASJL(J,LR,3)+SRHRS(I,J,LR)*COSZ 5422.
1174     740 ASJL(J,LR,4)=ASJL(J,LR,4)+TRHRS(I,J,LR) 5423.
1175     DJ(JR,1)=DJ(JR,1)+(S0*COSZ)*DXYPJ 5440.
1176     DJ(JR,2)=DJ(JR,2)+(SNFS(I,J,4)*COSZ)*DXYPJ 5444.
1177     DJ(JR,3)=DJ(JR,3)+(SNFS(I,J,1)*COSZ)*DXYPJ 5448.
1178     DJ(JR,5)=DJ(JR,5)+(SRHR(I,J,1)*COSZ/(ALB(I,J,1)+1.E-20))*DXYPJ 5452.
1179     DJ(JR,6)=DJ(JR,6)+(SRHR(I,J,1)*COSZ)*DXYPJ 5456.
1180     DJ(JR,55)=DJ(JR,55)+BTMPW(I,J)*DXYPJ 5460.
1181     DJ(JR,67)=DJ(JR,67)+TRINCG(I,J)*DXYPJ 5464.
1182     DJ(JR,70)=DJ(JR,70)-TNFS(I,J,4)*DXYPJ 5468.
1183     DJ(JR,71)=DJ(JR,71)-TNFS(I,J,1)*DXYPJ 5472.
1184     770 CONTINUE 5485.
1185     780 CONTINUE 5492.
1186     C**** 5504.
1187     C**** UPDATE THE TEMPERATURES BY RADIATION 5505.
1188     C**** 5506.
1189     800 DO 820 J=1,JM 5507.
1190     IMAX=IM 5508.
1191     IF(J.EQ.1.OR.J.EQ.JM) IMAX=1 5509.
1192     DO 820 LR=1,3 5510.
1193     DO 820 I=1,IMAX 5511.
1194     820 RQT(I,J,LR)=RQT(I,J,LR)+(SRHRS(I,J,LR)*COSZ2(I,J) 5512.
1195     * +TRHRS(I,J,LR))*COE(LR+LM) 5513.
1196     840 DO 860 J=1,JM 5514.
1197     #if ( defined HR_DATA )
1198     dswhr(j)=DSWSRF(j)*COSZ1(1,j)
1199     dlwhr(j)=DLWSRF(j)
1200     #endif
1201     #if ( defined CLM )
1202 jscott 1.3 i=1
1203     dsw4clm(i,j)=DSWSRF(j)*COSZ1(1,j)
1204     dlw4clm(i,j)=DLWSRF(j)
1205     swinr4clm(i,j)=DSWNIR(j)*COSZ1(1,j)
1206     swvis4clm(i,j)=DSWVIS(j)*COSZ1(1,j)
1207     ! if(TAU.gt.3633.0.and.TAU.lt.3744.0.and.J.eq.29) then
1208     if(J.eq.-29) then
1209     ! print *,TAU,j
1210     ! print *,'DLWSRF(j)=',DLWSRF(j),
1211     ! & 'dlw4clm(i,j)=',dlw4clm(i,j)
1212     write (229,*)TAU,TOFAY,dsw4clm(i,j),dlw4clm(i,j)
1213     endif
1214     if(J.eq.-30) then
1215     write (230,*)TAU,TOFAY,dsw4clm(i,j),dlw4clm(i,j)
1216     endif
1217 jscott 1.1 c For TEM
1218     swtd4tem(j)=swtd4tem(j)+S0*COSZ1(1,j)
1219     swsd4tem(j)=swsd4tem(j)+DSWSRF(j)*COSZ1(1,j)
1220     nradd4tem(j)=nradd4tem(j)+1
1221     #endif
1222 jscott 1.3 #if ( defined OCEAN_3D )
1223     solarinc_ocean(J)=solarinc_ocean(J)+SWIN(j,1)*COSZ1(1,j)
1224     solarnet_ocean(J)=solarnet_ocean(J)+SWNET(j,1)*COSZ1(1,j)
1225     solarinc_ice(J)=solarinc_ice(J)+SWIN(j,2)*COSZ1(1,j)
1226     solarnet_ice(J)=solarnet_ice(J)+SWNET(j,2)*COSZ1(1,j)
1227     navrado(j)=navrado(j)+1
1228     navrad(j)=navrad(j)+1
1229     #endif
1230 jscott 1.1 IMAX=IM 5515.
1231     IF(J.EQ.1.OR.J.EQ.JM) IMAX=1 5516.
1232     if(HPRNT)then
1233     if(J.eq.1)then
1234     print *,' radia TAU=',TAU
1235     print *,' Before 860'
1236     print *,(T(1,1,L),L=1,LM)
1237     endif
1238     endif
1239     DO 860 L=1,LM 5517.
1240     DO 860 I=1,IMAX 5518.
1241     860 T(I,J,L)=T(I,J,L)+(SRHR(I,J,L+1)*COSZ1(I,J)+TRHR(I,J,L+1)) 5519.
1242     * *COE(L)/(P(I,J)*PK(I,J,L)) 5520.
1243     if(HPRNT)then
1244     print *,' radia TAU=',TAU
1245     print *,' after 860'
1246     print 'B66',T(1,1,1)
1247     print *,'COSZ1(1,1)',COSZ1(1,1),' COE(1)=',COE(1)
1248     print *,'SRHR'
1249     print *,(SRHR(1,1,L+1),L=1,LM)
1250     print *,'TRHR'
1251     print *,(TRHR(1,1,L+1),L=1,LM)
1252     print *,'T'
1253     print *,(T(1,1,L),L=1,LM)
1254     endif
1255     c if(ncallclm.ge.9)stop
1256     RETURN 5521.
1257     END 5522.

  ViewVC Help
Powered by ViewVC 1.1.22