File indexing completed on 2018-03-02 18:40:18 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
8b150b4af9 Andr*0002 subroutine lwrio (nymd,nhms,bi,bj,myid,istrip,npcs,
0003 . low_level,mid_level,
082e38725b Andr*0004 . im,jm,lm,
5cede9ba3f Andr*0005 . pz,plz,plze,dpres,pkht,pkz,tz,qz,oz,
0006 . co2,cfc11,cfc12,cfc22,methane,n2o,emissivity,
1662f365b2 Andr*0007 . tgz,radlwg,st4,dst4,
0008 . dtradlw,dlwdtg,dtradlwc,lwgclr,
032fb71841 Andr*0009 . nlwcld,cldlw,clwmo,nlwlz,lwlz,
7d59649d1b Andr*0010 . lpnt,imstturb,qliqave,fccave,landtype)
1662f365b2 Andr*0011
0012 implicit none
0013
0014
0015
8b150b4af9 Andr*0016 integer nymd,nhms,istrip,npcs,bi,bj,myid
d2aaec7e02 Andr*0017 integer mid_level,low_level
1662f365b2 Andr*0018 integer im,jm,lm
a456aa407c Andr*0019 _RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1)
0020 _RL dpres(im,jm,lm),pkht(im,jm,lm+1),pkz(im,jm,lm)
0021 _RL tz(im,jm,lm),qz(im,jm,lm),oz(im,jm,lm)
0022 _RL co2,cfc11,cfc12,cfc22,methane(lm),n2o(lm)
0023 _RL emissivity(im,jm,10)
0024 _RL tgz(im,jm),radlwg(im,jm),st4(im,jm),dst4(im,jm)
0025 _RL dtradlw(im,jm,lm),dlwdtg (im,jm,lm)
0026 _RL dtradlwc(im,jm,lm),lwgclr(im,jm)
5cede9ba3f Andr*0027 integer nlwcld,nlwlz
a456aa407c Andr*0028 _RL cldlw(im,jm,lm),clwmo(im,jm,lm),lwlz(im,jm,lm)
5cede9ba3f Andr*0029 logical lpnt
0030 integer imstturb
a456aa407c Andr*0031 _RL qliqave(im,jm,lm),fccave(im,jm,lm)
5cede9ba3f Andr*0032 integer landtype(im,jm)
1662f365b2 Andr*0033
0034
0035
0036 integer i,j,l,n,nn
0037
a456aa407c Andr*0038 _RL cldtot (im,jm,lm)
0039 _RL cldmxo (im,jm,lm)
0040
0041 _RL pl(istrip,lm)
0042 _RL pk(istrip,lm)
0043 _RL pke(istrip,lm)
0044 _RL ple(istrip,lm+1)
0045
0046 _RL ADELPL(ISTRIP,lm)
0047 _RL dtrad(istrip,lm),dtradc(istrip,lm)
0048 _RL OZL(ISTRIP,lm),TZL(ISTRIP,lm)
0049 _RL SHZL(ISTRIP,lm),CLRO(ISTRIP,lm)
0050 _RL CLMO(ISTRIP,lm)
0051 _RL flx(ISTRIP,lm+1),flxclr(ISTRIP,lm+1)
0052 _RL cldlz(istrip,lm)
0053 _RL dfdts(istrip,lm+1),dtdtg(istrip,lm)
0054
0055 _RL emiss(istrip,10)
0056 _RL taual(istrip,lm,10)
0057 _RL ssaal(istrip,lm,10)
0058 _RL asyal(istrip,lm,10)
0059 _RL cwc(istrip,lm,3)
0060 _RL reff(istrip,lm,3)
0061 _RL tauc(istrip,lm,3)
0062
0063 _RL SGMT4(ISTRIP),TSURF(ISTRIP),dsgmt4(ISTRIP)
5cede9ba3f Andr*0064 integer lwi(istrip)
0065
032fb71841 Andr*0066 _RL tmpstrip(istrip,lm)
0067 _RL tmpimjm(im,jm,lm)
f1eddbb0f1 Andr*0068 _RL tempor1(im,jm),tempor2(im,jm)
032fb71841 Andr*0069
a456aa407c Andr*0070 _RL getcon,secday,convrt
8b150b4af9 Andr*0071 #ifdef ALLOW_DIAGNOSTICS
0072 logical diagnostics_is_on
0073 external diagnostics_is_on
0074 _RL tmpdiag(im,jm)
0075 #endif
5cede9ba3f Andr*0076
0077 logical high, trace, cldwater
032fb71841 Andr*0078
0079
0080 data high /.false./
0081 data trace /.false./
1662f365b2 Andr*0082 data cldwater /.false./
0083
0084
0085
0086
0087
0088 SECDAY = GETCON('SDAY')
0089 CONVRT = GETCON('GRAVITY') / ( 100.0 * GETCON('CP') )
0090
0091
0092
0093 if(imstturb.ne.0) then
0094 do L =1,lm
0095 do j =1,jm
0096 do i =1,im
e7070f537c Cons*0097 cldtot(i,j,L)=min(1.0 _d 0,max(cldlw(i,j,L),fccave(i,j,L)/
0098 $ imstturb))
0099 cldmxo(i,j,L) = min( 1.0 _d 0, clwmo(i,j,L) )
1662f365b2 Andr*0100 lwlz(i,j,L) = lwlz(i,j,L) + qliqave(i,j,L)/imstturb
0101 enddo
0102 enddo
0103 enddo
0104 else
0105 do L =1,lm
0106 do j =1,jm
0107 do i =1,im
32362b8fa8 Cons*0108 cldtot(i,j,L) = min( 1.0 _d 0,cldlw(i,j,L) )
0109 cldmxo(i,j,L) = min( 1.0 _d 0,clwmo(i,j,L) )
1662f365b2 Andr*0110 enddo
0111 enddo
0112 enddo
0113 endif
0114
0115
0116
0117
0118
0119 DO 1000 NN=1,NPCS
0120
0121
0122
0123
0124
0125 CALL STRIP (OZ,OZL ,im*jm,ISTRIP,lm ,NN)
0126 CALL STRIP (PLZE,ple ,im*jm,ISTRIP,lm+1 ,NN)
0127 CALL STRIP (PLZ ,pl ,im*jm,ISTRIP,lm ,NN)
0128 CALL STRIP (PKZ ,pk ,im*jm,ISTRIP,lm ,NN)
0129 CALL STRIP (PKHT,pke ,im*jm,ISTRIP,lm ,NN)
0130 CALL STRIP (TZ,TZL ,im*jm,ISTRIP,lm ,NN)
0131 CALL STRIP (qz,SHZL ,im*jm,ISTRIP,lm ,NN)
0132 CALL STRIP (cldtot,CLRO ,im*jm,ISTRIP,lm ,NN)
0133 CALL STRIP (cldmxo,CLMO ,im*jm,ISTRIP,lm ,NN)
0134 CALL STRIP (lwlz,cldlz ,im*jm,ISTRIP,lm ,NN)
0135 CALL STRIP (tgz,tsurf ,im*jm,ISTRIP,1 ,NN)
0136
0137 CALL STRIP (emissivity,emiss,im*jm,ISTRIP,10,NN)
0138
0139 call stripitint (landtype,lwi,im*jm,im*jm,istrip,1,nn)
0140
032fb71841 Andr*0141 DO L = 1,lm
0142 DO I = 1,ISTRIP
0143 ADELPL(I,L) = convrt / ( ple(I,L+1)-ple(I,L) )
0144 ENDDO
1662f365b2 Andr*0145 ENDDO
0146
0147
0148
0149 IF(NLWCLD.NE.0)THEN
0150 DO L = 1,lm
0151 DO I = 1,ISTRIP
32362b8fa8 Cons*0152 CLRO(I,L) = min( 1.0 _d 0,clro(i,L) )
0153 CLMO(I,L) = min( 1.0 _d 0,clmo(i,L) )
1662f365b2 Andr*0154 ENDDO
0155 ENDDO
0156 ENDIF
0157
124bb68fab Andr*0158
0159
1662f365b2 Andr*0160 DO L = 1,lm
0161 DO I = 1,ISTRIP
0162 TZL(I,L) = TZL(I,L)*pk(I,L)
0163 ENDDO
0164 ENDDO
0165 DO I = 1,ISTRIP
0166
0167 TZL(I,lm) = ( 2*tzl(i,lm)+tsurf(i) )/3.0
0168 ENDDO
0169
0170
0171
0172 call opthk ( tzl,pl,ple,cldlz,clro,clmo,lwi,istrip,1,lm,tauc )
0173
0174 do n = 1,3
0175 do L = 1,lm
0176 do i = 1,istrip
0177 tauc(i,L,n) = tauc(i,L,n)*0.75
0178 enddo
0179 enddo
0180 enddo
0181
0182
0183
0184
0185 do n = 1,10
0186 do L = 1,lm
0187 do i = 1,istrip
0188 taual(i,L,n) = 0.
0189 ssaal(i,L,n) = 0.
0190 asyal(i,L,n) = 0.
0191 enddo
0192 enddo
0193 enddo
0194
0195
0196
0197 do n = 1,3
0198 do L = 1,lm
0199 do i = 1,istrip
0200 cwc(i,L,n) = 0.
0201 reff(i,L,n) = 0.
0202 enddo
0203 enddo
0204 enddo
0205
0206
0207
0208
0209
0210 call irrad ( istrip,1,1,lm,ple,tzl,shzl,ozl,tsurf,co2,
0211 . n2o,methane,cfc11,cfc12,cfc22,emiss,
0212 . cldwater,cwc,tauc,reff,clro,mid_level,low_level,
0213 . taual,ssaal,asyal,
0214 . high,trace,flx,flxclr,dfdts,sgmt4 )
0215
0216
0217
0218
0219
0220 do L = 1,lm
0221 do i = 1,istrip
0222 dtrad(i,L) = ( flx(i,L)- flx(i,L+1))*adelpl(i,L)
154b1b9078 Andr*0223 tmpstrip(i,L) = flx(i,L)
1662f365b2 Andr*0224 dtdtg(i,L) = ( dfdts(i,L)- dfdts(i,L+1))*adelpl(i,L)
0225 dtradc(i,L) = (flxclr(i,L)-flxclr(i,L+1))*adelpl(i,L)
0226 enddo
0227 enddo
0228
0229
0230
0231
0232
0233 DO I = 1,ISTRIP
0234 flx (i,1) = -flx (i,1)
0235 flxclr(i,1) = -flxclr(i,1)
0236 flx (i,lm+1) = -flx (i,lm+1)
0237 flxclr(i,lm+1) = -flxclr(i,lm+1)
0238 sgmt4(i) = - sgmt4(i)
0239 dsgmt4(i) = - dfdts(i,lm+1)
0240 ENDDO
0241
0242 CALL PASTE ( flx (1,lm+1),RADLWG,ISTRIP,im*jm,1,NN )
0243 CALL PASTE ( flxclr(1,lm+1),LWGCLR,ISTRIP,im*jm,1,NN )
0244
0245 CALL PASTE ( sgmt4, st4,ISTRIP,im*jm,1,NN )
0246 CALL PASTE ( dsgmt4,dst4,ISTRIP,im*jm,1,NN )
0247
0248
0249
0250
0251
f1eddbb0f1 Andr*0252 CALL PASTE(flx(1,1),tempor1,ISTRIP,im*jm,1,NN)
0253 CALL PASTE(flxclr(1,1),tempor2,ISTRIP,im*jm,1,NN)
1662f365b2 Andr*0254
0255
0256
0257
0258
0259 DO L = 1,lm
0260 DO I = 1,ISTRIP
032fb71841 Andr*0261 DTRAD (I,L) = ple(i,lm+1) * DTRAD (I,L)/pk(I,L)
0262 DTRADC(I,L) = ple(i,lm+1) * DTRADC(I,L)/pk(I,L)
0263 dtdtg(I,L) = ple(i,lm+1) * dtdtg (I,L)/pk(I,L)
1662f365b2 Andr*0264 ENDDO
0265 ENDDO
032fb71841 Andr*0266 CALL PASTE ( tmpstrip ,tmpimjm ,ISTRIP,im*jm,lm,NN )
1662f365b2 Andr*0267 CALL PASTE ( DTRAD ,DTRADLW ,ISTRIP,im*jm,lm,NN )
0268 CALL PASTE ( DTRADC,DTRADLWC,ISTRIP,im*jm,lm,NN )
0269 CALL PASTE ( dtdtg ,dlwdtg ,ISTRIP,im*jm,lm,NN )
0270
0271 1000 CONTINUE
0272
0273
0274
0275
0276
8b150b4af9 Andr*0277 #ifdef ALLOW_DIAGNOSTICS
0278 if(diagnostics_is_on('TGRLW ',myid) ) then
0279 call diagnostics_fill(tgz,'TGRLW ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0280 endif
0281
0282 do L = 1,lm
0283
8b150b4af9 Andr*0284 if(diagnostics_is_on('TLW ',myid) ) then
0285 do j = 1,jm
0286 do i = 1,im
0287 tmpdiag(i,j) = tz(i,j,L)*pkz(i,j,L)
0288 enddo
0289 enddo
0290 call diagnostics_fill(tmpdiag,'TLW ',L,1,3,bi,bj,myid)
0291 endif
1662f365b2 Andr*0292
8b150b4af9 Andr*0293 if(diagnostics_is_on('SHRAD ',myid) ) then
0294 do j = 1,jm
0295 do i = 1,im
0296 tmpdiag(i,j) = qz(i,j,L)*1000.
0297 enddo
0298 enddo
0299 call diagnostics_fill(tmpdiag,'SHRAD ',L,1,3,bi,bj,myid)
0300 endif
0301
0302 if(diagnostics_is_on('OZLW ',myid) ) then
0303 do j = 1,jm
0304 do i = 1,im
0305 tmpdiag(i,j) = oz(i,j,L)
0306 enddo
0307 enddo
0308 call diagnostics_fill(tmpdiag,'OZLW ',L,1,3,bi,bj,myid)
0309 endif
032fb71841 Andr*0310
f1eddbb0f1 Andr*0311 enddo
8b150b4af9 Andr*0312
0313 if(diagnostics_is_on('OLR ',myid) ) then
0314 call diagnostics_fill(tempor1,'OLR ',0,1,3,bi,bj,myid)
f1eddbb0f1 Andr*0315 endif
0316
8b150b4af9 Andr*0317 if(diagnostics_is_on('OLRCLR ',myid) ) then
0318 call diagnostics_fill(tempor2,'OLRCLR ',0,1,3,bi,bj,myid)
032fb71841 Andr*0319 endif
8b150b4af9 Andr*0320 #endif
1662f365b2 Andr*0321
0322
0323
0324
0325
0326 nlwlz = 0
0327 nlwcld = 0
0328 imstturb = 0
0329
0330 do L = 1,lm
0331 do j = 1,jm
0332 do i = 1,im
0333 fccave(i,j,L) = 0.
0334 qliqave(i,j,L) = 0.
0335 enddo
0336 enddo
0337 enddo
0338
0339 return
0340 end
0341
0342
0343 subroutine irrad (m,n,ndim,np,pl,ta,wa,oa,ts,co2,
0344 * n2o,ch4,cfc11,cfc12,cfc22,emiss,
0345 * cldwater,cwc,taucl,reff,fcld,ict,icb,
0346 * taual,ssaal,asyal,
0347 * high,trace,flx,flc,dfdts,st4)
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446
0447
0448
0449
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468
0469
0470
0471
0472
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523
0524
0525
0526
0527
0528
0529
0530
0531
0532
0533
0534
0535
0536
0537
0538
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548
0549
0550
0551
0552
0553
0554
0555
0556
0557
0558 implicit none
0559 integer nx,no,nc,nh,nt
0560 integer i,j,k,ip,iw,it,ib,ik,iq,isb,k1,k2
0561 parameter (nx=26,no=21,nc=24,nh=31,nt=7)
0562
0563
0564
0565 integer m,n,ndim,np,ict,icb
a456aa407c Andr*0566 _RL pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np),
1662f365b2 Andr*0567 * ts(m,ndim)
a456aa407c Andr*0568 _RL co2,n2o(np),ch4(np),cfc11,cfc12,cfc22,emiss(m,ndim,10)
0569 _RL cwc(m,ndim,np,3),taucl(m,ndim,np,3),reff(m,ndim,np,3),
1662f365b2 Andr*0570 * fcld(m,ndim,np)
a456aa407c Andr*0571 _RL taual(m,ndim,np,10),ssaal(m,ndim,np,10),asyal(m,ndim,np,10)
1662f365b2 Andr*0572 logical cldwater,high,trace
0573
0574
0575
a456aa407c Andr*0576 _RL flx(m,ndim,np+1),flc(m,ndim,np+1),dfdts(m,ndim,np+1),
1662f365b2 Andr*0577 * st4(m,ndim)
0578
0579
0580
a456aa407c Andr*0581 _RL cb(5,10)
0582 _RL xkw(9),aw(9),bw(9),pm(9),fkw(6,9),gkw(6,3),xke(9)
0583 _RL aib(3,10),awb(4,10),aiw(4,10),aww(4,10),aig(4,10),awg(4,10)
1662f365b2 Andr*0584 integer ne(9),mw(9)
0585
0586
0587
a456aa407c Andr*0588 _RL pa(m,n,np),dt(m,n,np)
0589 _RL sh2o(m,n,np+1),swpre(m,n,np+1),swtem(m,n,np+1)
0590 _RL sco3(m,n,np+1),scopre(m,n,np+1),scotem(m,n,np+1)
0591 _RL dh2o(m,n,np),dcont(m,n,np),dco2(m,n,np),do3(m,n,np)
0592 _RL dn2o(m,n,np),dch4(m,n,np)
0593 _RL df11(m,n,np),df12(m,n,np),df22(m,n,np)
0594 _RL th2o(m,n,6),tcon(m,n,3),tco2(m,n,6,2)
0595 _RL tn2o(m,n,4),tch4(m,n,4),tcom(m,n,2)
0596 _RL tf11(m,n),tf12(m,n),tf22(m,n)
0597 _RL h2oexp(m,n,np,6),conexp(m,n,np,3),co2exp(m,n,np,6,2)
0598 _RL n2oexp(m,n,np,4),ch4exp(m,n,np,4),comexp(m,n,np,2)
0599 _RL f11exp(m,n,np), f12exp(m,n,np), f22exp(m,n,np)
0600 _RL clr(m,n,0:np+1),fclr(m,n)
0601 _RL blayer(m,n,0:np+1),dlayer(m,n,np+1),dbs(m,n)
0602 _RL clrlw(m,n),clrmd(m,n),clrhi(m,n)
0603 _RL cwp(m,n,np,3)
0604 _RL trant(m,n),tranal(m,n),transfc(m,n,np+1),trantcr(m,n,np+1)
0605 _RL flxu(m,n,np+1),flxd(m,n,np+1),flcu(m,n,np+1),flcd(m,n,np+1)
0606 _RL rflx(m,n,np+1),rflc(m,n,np+1)
1662f365b2 Andr*0607
0608 logical oznbnd,co2bnd,h2otbl,conbnd,n2obnd
0609 logical ch4bnd,combnd,f11bnd,f12bnd,f22bnd,b10bnd
0610
a456aa407c Andr*0611 _RL c1 (nx,nc,nt),c2 (nx,nc,nt),c3 (nx,nc,nt)
0612 _RL o1 (nx,no,nt),o2 (nx,no,nt),o3 (nx,no,nt)
0613 _RL h11(nx,nh,nt),h12(nx,nh,nt),h13(nx,nh,nt)
0614 _RL h21(nx,nh,nt),h22(nx,nh,nt),h23(nx,nh,nt)
0615 _RL h81(nx,nh,nt),h82(nx,nh,nt),h83(nx,nh,nt)
1662f365b2 Andr*0616
a456aa407c Andr*0617 _RL dp,xx,p1,dwe,dpe,a1,b1,fk1,a2,b2,fk2
0618 _RL w1,w2,w3,g1,g2,g3,ww,gg,ff,taux,reff1,reff2
1662f365b2 Andr*0619
0620
0621
0622
0623
0624 data cb/
0625 1 -2.6844e-1,-8.8994e-2, 1.5676e-3,-2.9349e-6, 2.2233e-9,
0626 2 3.7315e+1,-7.4758e-1, 4.6151e-3,-6.3260e-6, 3.5647e-9,
0627 3 3.7187e+1,-3.9085e-1,-6.1072e-4, 1.4534e-5,-1.6863e-8,
0628 4 -4.1928e+1, 1.0027e+0,-8.5789e-3, 2.9199e-5,-2.5654e-8,
0629 5 -4.9163e+1, 9.8457e-1,-7.0968e-3, 2.0478e-5,-1.5514e-8,
0630 6 -4.7107e+1, 8.9130e-1,-5.9735e-3, 1.5596e-5,-9.5911e-9,
0631 7 -5.4041e+1, 9.5332e-1,-5.7784e-3, 1.2555e-5,-2.9377e-9,
0632 8 -6.9233e+0,-1.5878e-1, 3.9160e-3,-2.4496e-5, 4.9301e-8,
0633 9 1.1483e+2,-2.2376e+0, 1.6394e-2,-5.3672e-5, 6.6456e-8,
0634 * 1.9668e+1,-3.1161e-1, 1.2886e-3, 3.6835e-7,-1.6212e-9/
0635
0636
0637
0638
0639
0640 data xkw / 29.55 , 4.167e-1, 1.328e-2, 5.250e-4,
0641 * 5.25e-4, 9.369e-3, 4.719e-2, 1.320e-0, 5.250e-4/
0642
0643
0644
0645
0646
0647 data xke / 0.00, 0.00, 27.40, 15.8,
0648 * 9.40, 7.75, 0.0, 0.0, 0.0/
0649
0650
0651
0652
0653 data mw /6,6,8,6,6,8,9,6,16/
0654
0655
0656
0657
0658 data aw/ 0.0021, 0.0140, 0.0167, 0.0302,
0659 * 0.0307, 0.0195, 0.0152, 0.0008, 0.0096/
0660 data bw/ -1.01e-5, 5.57e-5, 8.54e-5, 2.96e-4,
0661 * 2.86e-4, 1.108e-4, 7.608e-5, -3.52e-6, 1.64e-5/
0662
0663 data pm/ 1.0, 1.0, 1.0, 1.0, 1.0, 0.77, 0.5, 1.0, 1.0/
0664
0665
0666
0667
0668
0669 data fkw / 0.2747,0.2717,0.2752,0.1177,0.0352,0.0255,
0670 2 0.1521,0.3974,0.1778,0.1826,0.0374,0.0527,
0671 3 6*1.00,
0672 4 0.4654,0.2991,0.1343,0.0646,0.0226,0.0140,
0673 5 0.5543,0.2723,0.1131,0.0443,0.0160,0.0000,
0674 6 0.5955,0.2693,0.0953,0.0335,0.0064,0.0000,
0675 7 0.1958,0.3469,0.3147,0.1013,0.0365,0.0048,
0676 8 0.0740,0.1636,0.4174,0.1783,0.1101,0.0566,
0677 9 0.1437,0.2197,0.3185,0.2351,0.0647,0.0183/
0678
0679
0680
0681
0682
0683
0684 data gkw/ 0.1782,0.0593,0.0215,0.0068,0.0022,0.0000,
0685 2 0.0923,0.1675,0.0923,0.0187,0.0178,0.0000,
0686 3 0.0000,0.1083,0.1581,0.0455,0.0274,0.0041/
0687
0688
0689
0690
0691 data ne /0,0,3,1,1,1,0,0,0/
0692
0693
0694
0695
0696
0697 data aib / -0.44171, 0.62951, 0.06465,
0698 2 -0.13727, 0.61291, 0.28962,
0699 3 -0.01878, 1.67680, 0.79080,
0700 4 -0.01896, 1.06510, 0.69493,
0701 5 -0.04788, 0.88178, 0.54492,
0702 6 -0.02265, 1.57390, 0.76161,
0703 7 -0.01038, 2.15640, 0.89045,
0704 8 -0.00450, 2.51370, 0.95989,
0705 9 -0.00044, 3.15050, 1.03750,
0706 * -0.02956, 1.44680, 0.71283/
0707
0708
0709
0710
0711
0712 data awb / 0.08641, 0.01769, -1.5572e-3, 3.4896e-5,
0713 2 0.22027, 0.00997, -1.8719e-3, 5.3112e-5,
0714 3 0.39252, -0.02817, 7.2931e-4, -3.8151e-6,
0715 4 0.39975, -0.03426, 1.2884e-3, -1.7930e-5,
0716 5 0.34021, -0.02805, 1.0654e-3, -1.5443e-5,
0717 6 0.15587, 0.00371, -7.7705e-4, 2.0547e-5,
0718 7 0.05518, 0.04544, -4.2067e-3, 1.0184e-4,
0719 8 0.12724, 0.04751, -5.2037e-3, 1.3711e-4,
0720 9 0.30390, 0.01656, -3.5271e-3, 1.0828e-4,
0721 * 0.63617, -0.06287, 2.2350e-3, -2.3177e-5/
0722
0723
0724
0725
0726
0727 data aiw/ 0.17201, 1.2229e-2, -1.4837e-4, 5.8020e-7,
0728 2 0.81470, -2.7293e-3, 9.7816e-8, 5.7650e-8,
0729 3 0.54859, -4.8273e-4, 5.4353e-6, -1.5679e-8,
0730 4 0.39218, 4.1717e-3, - 4.8869e-5, 1.9144e-7,
0731 5 0.71773, -3.3640e-3, 1.9713e-5, -3.3189e-8,
0732 6 0.77345, -5.5228e-3, 4.8379e-5, -1.5151e-7,
0733 7 0.74975, -5.6604e-3, 5.6475e-5, -1.9664e-7,
0734 8 0.69011, -4.5348e-3, 4.9322e-5, -1.8255e-7,
0735 9 0.83963, -6.7253e-3, 6.1900e-5, -2.0862e-7,
0736 * 0.64860, -2.8692e-3, 2.7656e-5, -8.9680e-8/
0737
0738
0739
0740
0741
0742 data aww/ -7.8566e-2, 8.0875e-2, -4.3403e-3, 8.1341e-5,
0743 2 -1.3384e-2, 9.3134e-2, -6.0491e-3, 1.3059e-4,
0744 3 3.7096e-2, 7.3211e-2, -4.4211e-3, 9.2448e-5,
0745 4 -3.7600e-3, 9.3344e-2, -5.6561e-3, 1.1387e-4,
0746 5 0.40212, 7.8083e-2, -5.9583e-3, 1.2883e-4,
0747 6 0.57928, 5.9094e-2, -5.4425e-3, 1.2725e-4,
0748 7 0.68974, 4.2334e-2, -4.9469e-3, 1.2863e-4,
0749 8 0.80122, 9.4578e-3, -2.8508e-3, 9.0078e-5,
0750 9 1.02340, -2.6204e-2, 4.2552e-4, 3.2160e-6,
0751 * 0.05092, 7.5409e-2, -4.7305e-3, 1.0121e-4/
0752
0753
0754
0755
0756 data aig / 0.57867, 1.0135e-2, -1.1142e-4, 4.1537e-7,
0757 2 0.72259, 3.1149e-3, -1.9927e-5, 5.6024e-8,
0758 3 0.76109, 4.5449e-3, -4.6199e-5, 1.6446e-7,
0759 4 0.86934, 2.7474e-3, -3.1301e-5, 1.1959e-7,
0760 5 0.89103, 1.8513e-3, -1.6551e-5, 5.5193e-8,
0761 6 0.86325, 2.1408e-3, -1.6846e-5, 4.9473e-8,
0762 7 0.85064, 2.5028e-3, -2.0812e-5, 6.3427e-8,
0763 8 0.86945, 2.4615e-3, -2.3882e-5, 8.2431e-8,
0764 9 0.80122, 3.1906e-3, -2.4856e-5, 7.2411e-8,
0765 * 0.73290, 4.8034e-3, -4.4425e-5, 1.4839e-7/
0766
0767
0768
0769
0770 data awg / -0.51930, 0.20290, -1.1747e-2, 2.3868e-4,
0771 2 -0.22151, 0.19708, -1.2462e-2, 2.6646e-4,
0772 3 0.14157, 0.14705, -9.5802e-3, 2.0819e-4,
0773 4 0.41590, 0.10482, -6.9118e-3, 1.5115e-4,
0774 5 0.55338, 7.7016e-2, -5.2218e-3, 1.1587e-4,
0775 6 0.61384, 6.4402e-2, -4.6241e-3, 1.0746e-4,
0776 7 0.67891, 4.8698e-2, -3.7021e-3, 9.1966e-5,
0777 8 0.78169, 2.0803e-2, -1.4749e-3, 3.9362e-5,
0778 9 0.93218, -3.3425e-2, 2.9632e-3, -6.9362e-5,
0779 * 0.01649, 0.16561, -1.0723e-2, 2.3220e-4/
0780
0781
0782
0783
0784 logical first
0785 data first /.true./
0786
65b2c6ac12 Andr*0787 #include "h2o-tran3.h"
0788 #include "co2-tran3.h"
0789 #include "o3-tran3.h"
1662f365b2 Andr*0790
032fb71841 Andr*0791
0792
1662f365b2 Andr*0793
154b1b9078 Andr*0794 if (first) then
1662f365b2 Andr*0795
0796
0797
0798 if (high) then
0799
0800 do iw=1,nh
0801 do ip=1,nx
0802 h11(ip,iw,1)=1.0-h11(ip,iw,1)
0803 h21(ip,iw,1)=1.0-h21(ip,iw,1)
0804 h81(ip,iw,1)=1.0-h81(ip,iw,1)
0805 enddo
0806 enddo
0807
0808 do iw=1,nc
0809 do ip=1,nx
0810 c1(ip,iw,1)=1.0-c1(ip,iw,1)
0811 enddo
0812 enddo
0813
0814
0815
0816
0817
0818
0819 do it=2,nt
0820 do iw=1,nc
0821 do ip=1,nx
0822 c1 (ip,iw,it)= c1(ip,iw,1)
0823 c2 (ip,iw,it)= c2(ip,iw,1)
0824 c3 (ip,iw,it)= c3(ip,iw,1)
0825 enddo
0826 enddo
0827 do iw=1,nh
0828 do ip=1,nx
0829 h11(ip,iw,it)=h11(ip,iw,1)
0830 h12(ip,iw,it)=h12(ip,iw,1)
0831 h13(ip,iw,it)=h13(ip,iw,1)
0832 h21(ip,iw,it)=h21(ip,iw,1)
0833 h22(ip,iw,it)=h22(ip,iw,1)
0834 h23(ip,iw,it)=h23(ip,iw,1)
0835 h81(ip,iw,it)=h81(ip,iw,1)
0836 h82(ip,iw,it)=h82(ip,iw,1)
0837 h83(ip,iw,it)=h83(ip,iw,1)
0838 enddo
0839 enddo
0840 enddo
0841
0842 endif
0843
0844
0845
0846 do iw=1,no
0847 do ip=1,nx
0848 o1(ip,iw,1)=1.0-o1(ip,iw,1)
0849 enddo
0850 enddo
0851
0852 do it=2,nt
0853 do iw=1,no
0854 do ip=1,nx
0855 o1 (ip,iw,it)= o1(ip,iw,1)
0856 o2 (ip,iw,it)= o2(ip,iw,1)
0857 o3 (ip,iw,it)= o3(ip,iw,1)
0858 enddo
0859 enddo
0860 enddo
0861
154b1b9078 Andr*0862 first=.false.
1662f365b2 Andr*0863
154b1b9078 Andr*0864 endif
1662f365b2 Andr*0865
0866
0867
0868
0869 do j=1,n
0870 do i=1,m
0871 if (pl(i,j,1) .eq. 0.0) then
0872 pl(i,j,1)=1.0e-4
0873 endif
0874 enddo
0875 enddo
0876
0877
0878
0879 do k=1,np
0880 do j=1,n
0881 do i=1,m
0882 pa(i,j,k)=0.5*(pl(i,j,k)+pl(i,j,k+1))
0883 dt(i,j,k)=ta(i,j,k)-250.0
0884 enddo
0885 enddo
0886 enddo
0887
0888
0889
0890
0891
0892
0893
0894
0895
0896
0897
0898
0899
0900
0901
0902
0903
0904
0905
0906 do k=1,np
0907 do j=1,n
0908 do i=1,m
0909
0910 dp = pl(i,j,k+1)-pl(i,j,k)
0911 dh2o (i,j,k) = 1.02*wa(i,j,k)*dp+1.e-10
0912 do3 (i,j,k) = 476.*oa(i,j,k)*dp+1.e-10
0913 dco2 (i,j,k) = 789.*co2*dp+1.e-10
0914 dch4 (i,j,k) = 789.*ch4(k)*dp+1.e-10
0915 dn2o (i,j,k) = 789.*n2o(k)*dp+1.e-10
0916 df11 (i,j,k) = 789.*cfc11*dp+1.e-10
0917 df12 (i,j,k) = 789.*cfc12*dp+1.e-10
0918 df22 (i,j,k) = 789.*cfc22*dp+1.e-10
0919
0920
0921
0922
0923 xx=pa(i,j,k)*0.001618*wa(i,j,k)*wa(i,j,k)*dp
0924 dcont(i,j,k) = xx*exp(1800./ta(i,j,k)-6.081)+1.e-10
0925
0926 enddo
0927 enddo
0928 enddo
0929
0930
0931
0932
0933 if (high) then
0934 call column(m,n,np,pa,dt,dh2o,sh2o,swpre,swtem)
0935 endif
0936
0937
0938
0939
0940
0941
0942 if (cldwater) then
0943 do k=1,np
0944 do j=1,n
0945 do i=1,m
0946 dp =pl(i,j,k+1)-pl(i,j,k)
0947 cwp(i,j,k,1)=1.02*10000.0*cwc(i,j,k,1)*dp
0948 cwp(i,j,k,2)=1.02*10000.0*cwc(i,j,k,2)*dp
0949 cwp(i,j,k,3)=1.02*10000.0*cwc(i,j,k,3)*dp
0950 taucl(i,j,k,3)=0.00307*cwp(i,j,k,3)
0951 enddo
0952 enddo
0953 enddo
0954 endif
0955
0956
0957
0958
0959
0960
0961
0962 do j=1,n
0963 do i=1,m
0964 clr(i,j,0) = 1.0
0965 clr(i,j,np+1) = 0.0
0966 st4(i,j) = 0.0
0967 transfc(i,j,np+1)=1.0
0968 trantcr(i,j,np+1)=1.0
0969 enddo
0970 enddo
0971
0972
0973
0974 do k=1,np+1
0975 do j=1,n
0976 do i=1,m
0977 flx(i,j,k) = 0.0
0978 flc(i,j,k) = 0.0
0979 dfdts(i,j,k)= 0.0
0980 rflx(i,j,k) = 0.0
0981 rflc(i,j,k) = 0.0
0982 enddo
0983 enddo
0984 enddo
0985
0986
0987
0988 do 1000 ib=1,10
0989
0990
0991
0992
0993
0994
0995
0996
0997
0998
0999
1000
1001
1002 h2otbl=high.and.(ib.eq.1.or.ib.eq.2.or.ib.eq.8)
1003 conbnd=ib.ge.3.and.ib.le.6
1004 co2bnd=ib.eq.3
1005 oznbnd=ib.eq.5
1006 n2obnd=ib.eq.6.or.ib.eq.7
1007 ch4bnd=ib.eq.6.or.ib.eq.7
1008 combnd=ib.eq.4.or.ib.eq.5
1009 f11bnd=ib.eq.4.or.ib.eq.5
1010 f12bnd=ib.eq.4.or.ib.eq.6
1011 f22bnd=ib.eq.4.or.ib.eq.6
1012 b10bnd=ib.eq.10
1013
1014
1015
1016
1017
1018 do k=1,np
1019 do j=1,n
1020 do i=1,m
1021 blayer(i,j,k)=ta(i,j,k)*(ta(i,j,k)*(ta(i,j,k)
1022 * *(ta(i,j,k)*cb(5,ib)+cb(4,ib))+cb(3,ib))
1023 * +cb(2,ib))+cb(1,ib)
1024 enddo
1025 enddo
1026 enddo
1027
3daafce20b Jean*1028
1662f365b2 Andr*1029
1030 do j=1,n
1031 do i=1,m
1032 blayer(i,j,0) = 0.0
1033 blayer(i,j,np+1)= ( ts(i,j)*(ts(i,j)*(ts(i,j)
1034 * *(ts(i,j)*cb(5,ib)+cb(4,ib))+cb(3,ib))
1035 * +cb(2,ib))+cb(1,ib) )*emiss(i,j,ib)
1036
1037
1038
1039
1040 dbs(i,j)=(ts(i,j)*(ts(i,j)*(ts(i,j)*4.*cb(5,ib)
1041 * +3.*cb(4,ib))+2.*cb(3,ib))+cb(2,ib))*emiss(i,j,ib)
1042
1043 enddo
1044 enddo
1045
1046 do k=1,np+1
1047 do j=1,n
1048 do i=1,m
1049 dlayer(i,j,k)=blayer(i,j,k-1)-blayer(i,j,k)
1050 enddo
1051 enddo
1052 enddo
1053
1054
1055
1056
1057
1058
1059
1060 if (high .and. co2bnd) then
1061 call column(m,n,np,pa,dt,dco2,sco3,scopre,scotem)
1062 endif
1063
1064 if (oznbnd) then
1065 call column(m,n,np,pa,dt,do3,sco3,scopre,scotem)
1066 endif
1067
1068
1069
1070 if (cldwater) then
1071 do k= 1, np
1072 do j= 1, n
1073 do i= 1, m
1074 taucl(i,j,k,1)=cwp(i,j,k,1)*(aib(1,ib)+aib(2,ib)/
1075 * reff(i,j,k,1)**aib(3,ib))
1076 taucl(i,j,k,2)=cwp(i,j,k,2)*(awb(1,ib)+(awb(2,ib)+(awb(3,ib)
1077 * +awb(4,ib)*reff(i,j,k,2))*reff(i,j,k,2))*reff(i,j,k,2))
1078 enddo
1079 enddo
1080 enddo
1081 endif
1082
1083
1084
1085
1086
1087
1088 do k= 1, np
1089 do j= 1, n
1090 do i= 1, m
1091
1092 clr(i,j,k) = 1.0
1093 taux=taucl(i,j,k,1)+taucl(i,j,k,2)+taucl(i,j,k,3)
1094
1095 if (taux.gt.0.02 .and. fcld(i,j,k).gt.0.01) then
1096
32362b8fa8 Cons*1097 reff1=min(reff(i,j,k,1),130. _d 0)
1098 reff2=min(reff(i,j,k,2),20.0 _d 0)
1662f365b2 Andr*1099
1100 w1=taucl(i,j,k,1)*(aiw(1,ib)+(aiw(2,ib)+(aiw(3,ib)
1101 * +aiw(4,ib)*reff1)*reff1)*reff1)
1102 w2=taucl(i,j,k,2)*(aww(1,ib)+(aww(2,ib)+(aww(3,ib)
1103 * +aww(4,ib)*reff2)*reff2)*reff2)
1104 w3=taucl(i,j,k,3)*0.54
1105 ww=(w1+w2+w3)/taux
1106
1107 g1=w1*(aig(1,ib)+(aig(2,ib)+(aig(3,ib)
1108 * +aig(4,ib)*reff1)*reff1)*reff1)
1109 g2=w2*(awg(1,ib)+(awg(2,ib)+(awg(3,ib)
1110 * +awg(4,ib)*reff2)*reff2)*reff2)
1111 g3=w3*0.95
1112
1113 gg=(g1+g2+g3)/(w1+w2+w3)
1114
1115
1116
1117
1118 ff=0.5+(0.3739+(0.0076+0.1185*gg)*gg)*gg
1119 taux=taux*(1.-ww*ff)
1120
1121
1122
1123
1124
1125 clr(i,j,k)=1.0-(fcld(i,j,k)*(1.0-exp(-1.66*taux)))
1126
1127 endif
1128
1129 enddo
1130 enddo
1131 enddo
1132
1133
1134
1135
1136 if (.not.h2otbl .and. .not.b10bnd) then
1137 call h2oexps(ib,m,n,np,dh2o,pa,dt,xkw,aw,bw,pm,mw,h2oexp)
1138 endif
1139
1140
1141
1142
1143 if (conbnd) then
1144 call conexps(ib,m,n,np,dcont,xke,conexp)
1145 endif
1146
1147
1148
1149
1150 if (.not.high .and. co2bnd) then
1151 call co2exps(m,n,np,dco2,pa,dt,co2exp)
1152 endif
1153
1154
1155
1156 if (trace) then
1157
1158
1159
1160 if (n2obnd) then
1161 call n2oexps(ib,m,n,np,dn2o,pa,dt,n2oexp)
1162 endif
1163
1164
1165
1166 if (ch4bnd) then
1167 call ch4exps(ib,m,n,np,dch4,pa,dt,ch4exp)
1168 endif
1169
1170
1171
1172 if (combnd) then
1173 call comexps(ib,m,n,np,dco2,dt,comexp)
1174 endif
1175
1176
1177
1178 if (f11bnd) then
1179 a1 = 1.26610e-3
1180 b1 = 3.55940e-6
1181 fk1 = 1.89736e+1
1182 a2 = 8.19370e-4
1183 b2 = 4.67810e-6
1184 fk2 = 1.01487e+1
1185 call cfcexps(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df11,dt,f11exp)
1186 endif
1187
1188
1189
1190 if (f12bnd) then
1191 a1 = 8.77370e-4
1192 b1 =-5.88440e-6
1193 fk1 = 1.58104e+1
1194 a2 = 8.62000e-4
1195 b2 =-4.22500e-6
1196 fk2 = 3.70107e+1
1197 call cfcexps(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df12,dt,f12exp)
1198 endif
1199
1200
1201
1202 if (f22bnd) then
1203 a1 = 9.65130e-4
1204 b1 = 1.31280e-5
1205 fk1 = 6.18536e+0
1206 a2 =-3.00010e-5
1207 b2 = 5.25010e-7
1208 fk2 = 3.27912e+1
1209 call cfcexps(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df22,dt,f22exp)
1210 endif
1211
1212
1213
1214
1215 if (b10bnd) then
1216 call b10exps(m,n,np,dh2o,dcont,dco2,dn2o,pa,dt
1217 * ,h2oexp,conexp,co2exp,n2oexp)
1218 endif
1219
1220 endif
1221
1222
1223
1224
1225
1226
1227
1228 do k=1,np+1
1229 do j=1,n
1230 do i=1,m
1231 flxu(i,j,k) = 0.0
1232 flxd(i,j,k) = 0.0
1233 flcu(i,j,k) = 0.0
1234 flcd(i,j,k) = 0.0
1235 enddo
1236 enddo
1237 enddo
1238
1239 do 2000 k1=1,np
1240
1241
1242
1243
1244
1245
1246
1247 do j=1,n
1248 do i=1,m
1249 fclr(i,j) = 1.0
1250 clrlw(i,j) = 1.0
1251 clrmd(i,j) = 1.0
1252 clrhi(i,j) = 1.0
1253 tranal(i,j)= 1.0
1254 enddo
1255 enddo
1256
1257
1258
1259 if (.not. h2otbl) then
1260 do ik=1,6
1261 do j=1,n
1262 do i=1,m
1263 th2o(i,j,ik)=1.0
1264 enddo
1265 enddo
1266 enddo
1267 endif
1268
1269
1270
1271 if (conbnd) then
1272 do iq=1,3
1273 do j=1,n
1274 do i=1,m
1275 tcon(i,j,iq)=1.0
1276 enddo
1277 enddo
1278 enddo
1279 endif
1280
1281
1282
1283
1284
1285 if (.not.high .and. co2bnd) then
1286 do isb=1,2
1287 do ik=1,6
1288 do j=1,n
1289 do i=1,m
1290 tco2(i,j,ik,isb)=1.0
1291 enddo
1292 enddo
1293 enddo
1294 enddo
1295 endif
1296
1297
1298
1299 if (trace) then
1300
1301
1302
1303 if (n2obnd) then
1304 do ik=1,4
1305 do j=1,n
1306 do i=1,m
1307 tn2o(i,j,ik)=1.0
1308 enddo
1309 enddo
1310 enddo
1311 endif
1312
1313
1314
1315 if (ch4bnd) then
1316 do ik=1,4
1317 do j=1,n
1318 do i=1,m
1319 tch4(i,j,ik)=1.0
1320 enddo
1321 enddo
1322 enddo
1323 endif
1324
1325
1326
1327 if (combnd) then
1328 do ik=1,2
1329 do j=1,n
1330 do i=1,m
1331 tcom(i,j,ik)=1.0
1332 enddo
1333 enddo
1334 enddo
1335 endif
1336
1337
1338
1339 if (f11bnd) then
1340 do j=1,n
1341 do i=1,m
1342 tf11(i,j)=1.0
1343 enddo
1344 enddo
1345 endif
1346
1347
1348
1349 if (f12bnd) then
1350 do j=1,n
1351 do i=1,m
1352 tf12(i,j)=1.0
1353 enddo
1354 enddo
1355 endif
1356
1357
1358
1359 if (f22bnd) then
1360 do j=1,n
1361 do i=1,m
1362 tf22(i,j)=1.0
1363 enddo
1364 enddo
1365 endif
1366
1367
1368
1369 if (b10bnd) then
1370 do ik=1,6
1371 do j=1,n
1372 do i=1,m
1373 th2o(i,j,ik)=1.0
1374 tco2(i,j,ik,1)=1.0
1375 enddo
1376 enddo
1377 enddo
1378 do iq=1,3
1379 do j=1,n
1380 do i=1,m
1381 tcon(i,j,iq)=1.0
1382 enddo
1383 enddo
1384 enddo
1385 do ik=1,4
1386 do j=1,n
1387 do i=1,m
1388 tn2o(i,j,ik)=1.0
1389 enddo
1390 enddo
1391 enddo
1392 endif
1393
1394 endif
1395
1396
1397
1398 do 3000 k2=k1+1,np+1
1399
1400
1401
1402 do j=1,n
1403 do i=1,m
1404 trant(i,j)=1.0
1405 enddo
1406 enddo
1407
1408 if (h2otbl) then
1409 w1=-8.0
1410 p1=-2.0
1411 dwe=0.3
1412 dpe=0.2
1413
1414
1415
1416 if (ib.eq.1) then
1417 call tablup(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,
1418 * w1,p1,dwe,dpe,h11,h12,h13,trant)
1419 endif
1420 if (ib.eq.2) then
1421 call tablup(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,
1422 * w1,p1,dwe,dpe,h21,h22,h23,trant)
1423
1424 endif
1425 if (ib.eq.8) then
1426 call tablup(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,
1427 * w1,p1,dwe,dpe,h81,h82,h83,trant)
1428 endif
1429
1430 else
1431
1432
1433
1434 if (.not.b10bnd) then
1435 call h2okdis(ib,m,n,np,k2-1,fkw,gkw,ne,h2oexp,conexp,
1436 * th2o,tcon,trant)
1437 endif
1438
1439 endif
1440
1441 if (co2bnd) then
1442
1443 if (high) then
1444
1445
1446
1447 w1=-4.0
1448 p1=-2.0
1449 dwe=0.3
1450 dpe=0.2
1451 call tablup(k1,k2,m,n,np,nx,nc,nt,sco3,scopre,scotem,
1452 * w1,p1,dwe,dpe,c1,c2,c3,trant)
046c3794c1 Andr*1453
1662f365b2 Andr*1454 else
1455
1456
1457 call co2kdis(m,n,np,k2-1,co2exp,tco2,trant)
1458
1459 endif
1460
1461 endif
1462
1463
1464
1465 if (oznbnd) then
1466 w1=-6.0
1467 p1=-2.0
1468 dwe=0.3
1469 dpe=0.2
1470 call tablup(k1,k2,m,n,np,nx,no,nt,sco3,scopre,scotem,
1471 * w1,p1,dwe,dpe,o1,o2,o3,trant)
154b1b9078 Andr*1472
1662f365b2 Andr*1473 endif
1474
1475
1476
1477 if (trace) then
1478
1479
1480
1481 if (n2obnd) then
1482 call n2okdis(ib,m,n,np,k2-1,n2oexp,tn2o,trant)
1483 endif
1484
1485
1486
1487 if (ch4bnd) then
1488 call ch4kdis(ib,m,n,np,k2-1,ch4exp,tch4,trant)
1489 endif
1490
1491
1492
1493 if (combnd) then
1494 call comkdis(ib,m,n,np,k2-1,comexp,tcom,trant)
1495 endif
1496
1497
1498
1499 if (f11bnd) then
1500 call cfckdis(m,n,np,k2-1,f11exp,tf11,trant)
1501 endif
1502
1503
1504
1505 if (f12bnd) then
1506 call cfckdis(m,n,np,k2-1,f12exp,tf12,trant)
1507 endif
1508
1509
1510
1511 if (f22bnd) then
1512 call cfckdis(m,n,np,k2-1,f22exp,tf22,trant)
1513 endif
1514
1515
1516
1517
1518 if (b10bnd) then
1519 call b10kdis(m,n,np,k2-1,h2oexp,conexp,co2exp,n2oexp
1520 * ,th2o,tcon,tco2,tn2o,trant)
1521 endif
1522
1523 endif
1524
1525
1526
1527 do j=1,n
1528 do i=1,m
1529 ff=0.5+(0.3739+(0.0076+0.1185*asyal(i,j,k2-1,ib))
1530 * *asyal(i,j,k2-1,ib))*asyal(i,j,k2-1,ib)
1531 taux=taual(i,j,k2-1,ib)*(1.-ssaal(i,j,k2-1,ib)*ff)
1532 tranal(i,j)=tranal(i,j)*exp(-1.66*taux)
1533 trant (i,j)=trant(i,j) *tranal(i,j)
1534 enddo
1535 enddo
1536
1537
1538
1539
1540
1541
1542 do j=1,n
1543 do i=1,m
1544 if( k2 .le. ict ) then
1545 clrhi(i,j)=min(clr(i,j,k2-1),clrhi(i,j))
1546 elseif( k2 .gt. ict .and. k2 .le. icb ) then
1547 clrmd(i,j)=min(clr(i,j,k2-1),clrmd(i,j))
1548 elseif( k2 .gt. icb ) then
1549 clrlw(i,j)=min(clr(i,j,k2-1),clrlw(i,j))
1550 endif
1551 fclr(i,j)=clrlw(i,j)*clrmd(i,j)*clrhi(i,j)
1552
1553 enddo
1554 enddo
1555
1556
1557
1558
1559
1560
1561 if (.not. b10bnd) then
1562
1563 if (k2 .eq. k1+1) then
1564
1565 do j=1,n
1566 do i=1,m
1567
1568
1569
1570 flcu(i,j,k1)=flcu(i,j,k1)-blayer(i,j,k1)
1571 flcd(i,j,k2)=flcd(i,j,k2)+blayer(i,j,k1)
1572
1573
1574
1575 flxu(i,j,k1)=flxu(i,j,k1)-blayer(i,j,k1)
1576 flxd(i,j,k2)=flxd(i,j,k2)+blayer(i,j,k1)
1577
1578 enddo
1579 enddo
1580
1581 endif
1582
1583
1584
1585
1586 do j=1,n
1587 do i=1,m
1588 xx=trant(i,j)*dlayer(i,j,k2)
1589 flcu(i,j,k1) =flcu(i,j,k1)+xx
1590 flxu(i,j,k1) =flxu(i,j,k1)+xx*fclr(i,j)
1591 xx=trant(i,j)*dlayer(i,j,k1)
1592 flcd(i,j,k2) =flcd(i,j,k2)+xx
1593 flxd(i,j,k2) =flxd(i,j,k2)+xx*fclr(i,j)
1594 enddo
1595 enddo
1596
1597 else
1598
1599
1600
1601 if (trace) then
1602
1603 do j=1,n
1604 do i=1,m
1605 rflx(i,j,k1) = rflx(i,j,k1)+trant(i,j)*fclr(i,j)
1606 * *dlayer(i,j,k2)
1607 rflx(i,j,k2) = rflx(i,j,k2)+trant(i,j)*fclr(i,j)
1608 * *dlayer(i,j,k1)
1609 rflc(i,j,k1) = rflc(i,j,k1)+trant(i,j)
1610 * *dlayer(i,j,k2)
1611 rflc(i,j,k2) = rflc(i,j,k2)+trant(i,j)
1612 * *dlayer(i,j,k1)
1613 enddo
1614 enddo
1615
1616 endif
1617
1618 endif
1619
1620 3000 continue
1621
1622
1623
1624
1625 do j=1,n
1626 do i=1,m
1627 trantcr(i,j,k1) =trant(i,j)
1628 transfc(i,j,k1) =trant(i,j)*fclr(i,j)
1629 enddo
1630 enddo
1631
1632
1633
1634
1635 if (trace .or. (.not. b10bnd) ) then
1636
1637 do j=1,n
1638 do i=1,m
1639 dfdts(i,j,k1) =dfdts(i,j,k1)-dbs(i,j)*transfc(i,j,k1)
1640 enddo
1641 enddo
1642
1643 endif
1644
1645 2000 continue
1646
1647 if (.not. b10bnd) then
1648
1649
1650
1651 do j=1,n
1652 do i=1,m
1653 flcu(i,j,np+1)=flcu(i,j,np+1)-blayer(i,j,np+1)
1654 flxu(i,j,np+1)=flxu(i,j,np+1)-blayer(i,j,np+1)
1655 st4(i,j)=st4(i,j)-blayer(i,j,np+1)
1656 dfdts(i,j,np+1)=dfdts(i,j,np+1)-dbs(i,j)
1657 enddo
1658 enddo
1659
1660
1661
1662
1663 do k=1, np+1
1664 do j=1,n
1665 do i=1,m
1666 flcu(i,j,k)=flcu(i,j,k)-
1667 * flcd(i,j,np+1)*trantcr(i,j,k)*(1.-emiss(i,j,ib))
1668 flxu(i,j,k)=flxu(i,j,k)-
1669 * flxd(i,j,np+1)*transfc(i,j,k)*(1.-emiss(i,j,ib))
1670 enddo
1671 enddo
1672 enddo
1673
1674 endif
1675
1676
1677
1678 do k=1,np+1
1679 do j=1,n
1680 do i=1,m
1681 flc(i,j,k)=flc(i,j,k)+flcd(i,j,k)+flcu(i,j,k)
1682 flx(i,j,k)=flx(i,j,k)+flxd(i,j,k)+flxu(i,j,k)
1683 enddo
1684 enddo
1685 enddo
1686
1687 1000 continue
1688
1689
1690
1691 do k=1,np+1
1692 do j=1,n
1693 do i=1,m
1694 flc(i,j,k)=flc(i,j,k)+rflc(i,j,k)
1695 flx(i,j,k)=flx(i,j,k)+rflx(i,j,k)
1696 enddo
1697 enddo
1698 enddo
1699
1700 return
1701 end
1702
1703 subroutine column (m,n,np,pa,dt,sabs0,sabs,spre,stem)
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727 implicit none
1728 integer m,n,np,i,j,k
1729
1730
1731
a456aa407c Andr*1732 _RL pa(m,n,np),dt(m,n,np),sabs0(m,n,np)
1662f365b2 Andr*1733
1734
1735
a456aa407c Andr*1736 _RL sabs(m,n,np+1),spre(m,n,np+1),stem(m,n,np+1)
1662f365b2 Andr*1737
1738
1739 do j=1,n
1740 do i=1,m
1741 sabs(i,j,1)=0.0
1742 spre(i,j,1)=0.0
1743 stem(i,j,1)=0.0
1744 enddo
1745 enddo
1746
1747 do k=1,np
1748 do j=1,n
1749 do i=1,m
1750 sabs(i,j,k+1)=sabs(i,j,k)+sabs0(i,j,k)
1751 spre(i,j,k+1)=spre(i,j,k)+pa(i,j,k)*sabs0(i,j,k)
1752 stem(i,j,k+1)=stem(i,j,k)+dt(i,j,k)*sabs0(i,j,k)
1753 enddo
1754 enddo
1755 enddo
1756
1757 return
1758 end
1759
1760 subroutine h2oexps(ib,m,n,np,dh2o,pa,dt,xkw,aw,bw,pm,mw,h2oexp)
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782 implicit none
1783 integer ib,m,n,np,i,j,k,ik
1784
1785
1786
a456aa407c Andr*1787 _RL dh2o(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*1788
1789
1790
a456aa407c Andr*1791 _RL h2oexp(m,n,np,6)
1662f365b2 Andr*1792
1793
1794
1795 integer mw(9)
a456aa407c Andr*1796 _RL xkw(9),aw(9),bw(9),pm(9)
1662f365b2 Andr*1797
1798
1799
a456aa407c Andr*1800 _RL xh,xh1
1662f365b2 Andr*1801
1802
1803
1804
1805
1806
1807 do k=1,np
1808 do j=1,n
1809 do i=1,m
1810
1811
1812
1813
1814 xh = dh2o(i,j,k)*(pa(i,j,k)/500.)**pm(ib)
1815 1 * ( 1.+(aw(ib)+bw(ib)* dt(i,j,k))*dt(i,j,k) )
1816
1817
1818
1819
1820 h2oexp(i,j,k,1) = exp(-xh*xkw(ib))
1821
1822 enddo
1823 enddo
1824 enddo
1825
1826 do ik=2,6
1827
1828 if (mw(ib).eq.6) then
1829
1830 do k=1,np
1831 do j=1,n
1832 do i=1,m
1833 xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
1834 h2oexp(i,j,k,ik) = xh*xh*xh
1835 enddo
1836 enddo
1837 enddo
1838
1839 elseif (mw(ib).eq.8) then
1840
1841 do k=1,np
1842 do j=1,n
1843 do i=1,m
1844 xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
1845 xh = xh*xh
1846 h2oexp(i,j,k,ik) = xh*xh
1847 enddo
1848 enddo
1849 enddo
1850
1851 elseif (mw(ib).eq.9) then
1852
1853 do k=1,np
1854 do j=1,n
1855 do i=1,m
1856 xh=h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
1857 xh1 = xh*xh
1858 h2oexp(i,j,k,ik) = xh*xh1
1859 enddo
1860 enddo
1861 enddo
1862
1863 else
1864
1865 do k=1,np
1866 do j=1,n
1867 do i=1,m
1868 xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
1869 xh = xh*xh
1870 xh = xh*xh
1871 h2oexp(i,j,k,ik) = xh*xh
1872 enddo
1873 enddo
1874 enddo
1875
1876 endif
1877 enddo
1878
1879 return
1880 end
1881
1882 subroutine conexps(ib,m,n,np,dcont,xke,conexp)
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898 implicit none
1899 integer ib,m,n,np,i,j,k,iq
1900
1901
1902
a456aa407c Andr*1903 _RL dcont(m,n,np)
1662f365b2 Andr*1904
1905
1906
a456aa407c Andr*1907 _RL conexp(m,n,np,3)
1662f365b2 Andr*1908
1909
1910
a456aa407c Andr*1911 _RL xke(9)
1662f365b2 Andr*1912
1913
1914
1915 do k=1,np
1916 do j=1,n
1917 do i=1,m
1918 conexp(i,j,k,1) = exp(-dcont(i,j,k)*xke(ib))
1919 enddo
1920 enddo
1921 enddo
1922
1923 if (ib .eq. 3) then
1924
1925
1926
1927
1928
1929 do iq=2,3
1930 do k=1,np
1931 do j=1,n
1932 do i=1,m
1933 conexp(i,j,k,iq) = conexp(i,j,k,iq-1) *conexp(i,j,k,iq-1)
1934 enddo
1935 enddo
1936 enddo
1937 enddo
1938
1939 endif
1940
1941 return
1942 end
1943
1944 subroutine co2exps(m,n,np,dco2,pa,dt,co2exp)
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959 implicit none
1960 integer m,n,np,i,j,k
1961
1962
1963
a456aa407c Andr*1964 _RL dco2(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*1965
1966
1967
a456aa407c Andr*1968 _RL co2exp(m,n,np,6,2)
1662f365b2 Andr*1969
1970
1971
a456aa407c Andr*1972 _RL xc
1662f365b2 Andr*1973
1974
1975
1976 do k=1,np
1977 do j=1,n
1978 do i=1,m
1979
1980
1981
1982
1983 xc = dco2(i,j,k)*(pa(i,j,k)/300.0)**0.5
1984 1 *(1.+(0.0182+1.07e-4*dt(i,j,k))*dt(i,j,k))
1985
1986
1987
1988 co2exp(i,j,k,1,1)=exp(-xc*2.656e-5)
1989
1990 xc=co2exp(i,j,k,1,1)*co2exp(i,j,k,1,1)
1991 xc=xc*xc
1992 co2exp(i,j,k,2,1)=xc*xc
1993
1994 xc=co2exp(i,j,k,2,1)*co2exp(i,j,k,2,1)
1995 xc=xc*xc
1996 co2exp(i,j,k,3,1)=xc*xc
1997
1998 xc=co2exp(i,j,k,3,1)*co2exp(i,j,k,3,1)
1999 xc=xc*xc
2000 co2exp(i,j,k,4,1)=xc*xc
2001
2002 xc=co2exp(i,j,k,4,1)*co2exp(i,j,k,4,1)
2003 xc=xc*xc
2004 co2exp(i,j,k,5,1)=xc*xc
2005
2006 xc=co2exp(i,j,k,5,1)*co2exp(i,j,k,5,1)
2007 xc=xc*xc
2008 co2exp(i,j,k,6,1)=xc*xc
2009
2010
2011
2012
2013 xc = dco2(i,j,k)*(pa(i,j,k)/30.0)**0.85
2014 1 *(1.+(0.0042+2.00e-5*dt(i,j,k))*dt(i,j,k))
2015
2016 co2exp(i,j,k,1,2)=exp(-xc*2.656e-3)
2017
2018 xc=co2exp(i,j,k,1,2)*co2exp(i,j,k,1,2)
2019 xc=xc*xc
2020 co2exp(i,j,k,2,2)=xc*xc
2021
2022 xc=co2exp(i,j,k,2,2)*co2exp(i,j,k,2,2)
2023 xc=xc*xc
2024 co2exp(i,j,k,3,2)=xc*xc
2025
2026 xc=co2exp(i,j,k,3,2)*co2exp(i,j,k,3,2)
2027 xc=xc*xc
2028 co2exp(i,j,k,4,2)=xc*xc
2029
2030 xc=co2exp(i,j,k,4,2)*co2exp(i,j,k,4,2)
2031 xc=xc*xc
2032 co2exp(i,j,k,5,2)=xc*xc
2033
2034 xc=co2exp(i,j,k,5,2)*co2exp(i,j,k,5,2)
2035 xc=xc*xc
2036 co2exp(i,j,k,6,2)=xc*xc
2037
2038 enddo
2039 enddo
2040 enddo
2041
2042 return
2043 end
2044
2045 subroutine n2oexps(ib,m,n,np,dn2o,pa,dt,n2oexp)
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061 implicit none
2062 integer ib,m,n,np,i,j,k
2063
2064
2065
a456aa407c Andr*2066 _RL dn2o(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*2067
2068
2069
a456aa407c Andr*2070 _RL n2oexp(m,n,np,4)
1662f365b2 Andr*2071
2072
2073
a456aa407c Andr*2074 _RL xc,xc1,xc2
1662f365b2 Andr*2075
2076
2077
2078 do k=1,np
2079 do j=1,n
2080 do i=1,m
2081
2082
2083
2084 if (ib.eq.6) then
2085
2086 xc=dn2o(i,j,k)*(1.+(1.9297e-3+4.3750e-6*dt(i,j,k))*dt(i,j,k))
2087 n2oexp(i,j,k,1)=exp(-xc*6.31582e-2)
2088
2089 xc=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)*n2oexp(i,j,k,1)
2090 xc1=xc*xc
2091 xc2=xc1*xc1
2092 n2oexp(i,j,k,2)=xc*xc1*xc2
2093
2094
2095
2096 else
2097
2098 xc=dn2o(i,j,k)*(pa(i,j,k)/500.0)**0.48
2099 * *(1.+(1.3804e-3+7.4838e-6*dt(i,j,k))*dt(i,j,k))
2100 n2oexp(i,j,k,1)=exp(-xc*5.35779e-2)
2101
2102 xc=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)
2103 xc=xc*xc
2104 n2oexp(i,j,k,2)=xc*xc
2105 xc=n2oexp(i,j,k,2)*n2oexp(i,j,k,2)
2106 xc=xc*xc
2107 n2oexp(i,j,k,3)=xc*xc
2108 xc=n2oexp(i,j,k,3)*n2oexp(i,j,k,3)
2109 xc=xc*xc
2110 n2oexp(i,j,k,4)=xc*xc
2111
2112 endif
2113
2114 enddo
2115 enddo
2116 enddo
2117
2118 return
2119 end
2120
2121 subroutine ch4exps(ib,m,n,np,dch4,pa,dt,ch4exp)
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137 implicit none
2138 integer ib,m,n,np,i,j,k
2139
2140
2141
a456aa407c Andr*2142 _RL dch4(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*2143
2144
2145
a456aa407c Andr*2146 _RL ch4exp(m,n,np,4)
1662f365b2 Andr*2147
2148
2149
a456aa407c Andr*2150 _RL xc
1662f365b2 Andr*2151
2152
2153
2154 do k=1,np
2155 do j=1,n
2156 do i=1,m
2157
2158
2159
2160 if (ib.eq.6) then
2161
2162 xc=dch4(i,j,k)*(1.+(1.7007e-2+1.5826e-4*dt(i,j,k))*dt(i,j,k))
2163 ch4exp(i,j,k,1)=exp(-xc*5.80708e-3)
2164
2165
2166
2167 else
2168
2169 xc=dch4(i,j,k)*(pa(i,j,k)/500.0)**0.65
2170 * *(1.+(5.9590e-4-2.2931e-6*dt(i,j,k))*dt(i,j,k))
2171 ch4exp(i,j,k,1)=exp(-xc*6.29247e-2)
2172
2173 xc=ch4exp(i,j,k,1)*ch4exp(i,j,k,1)*ch4exp(i,j,k,1)
2174 xc=xc*xc
2175 ch4exp(i,j,k,2)=xc*xc
2176
2177 xc=ch4exp(i,j,k,2)*ch4exp(i,j,k,2)*ch4exp(i,j,k,2)
2178 xc=xc*xc
2179 ch4exp(i,j,k,3)=xc*xc
2180
2181 xc=ch4exp(i,j,k,3)*ch4exp(i,j,k,3)*ch4exp(i,j,k,3)
2182 xc=xc*xc
2183 ch4exp(i,j,k,4)=xc*xc
2184
2185 endif
2186
2187 enddo
2188 enddo
2189 enddo
2190
2191 return
2192 end
2193
2194 subroutine comexps(ib,m,n,np,dcom,dt,comexp)
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209 implicit none
2210 integer ib,m,n,np,i,j,k
2211
2212
2213
a456aa407c Andr*2214 _RL dcom(m,n,np),dt(m,n,np)
1662f365b2 Andr*2215
2216
2217
a456aa407c Andr*2218 _RL comexp(m,n,np,2)
1662f365b2 Andr*2219
2220
2221
a456aa407c Andr*2222 _RL xc,xc1,xc2
1662f365b2 Andr*2223
2224
2225
2226 do k=1,np
2227 do j=1,n
2228 do i=1,m
2229
2230
2231
2232 if (ib.eq.4) then
2233
2234 xc=dcom(i,j,k)*(1.+(3.5775e-2+4.0447e-4*dt(i,j,k))*dt(i,j,k))
2235 comexp(i,j,k,1)=exp(-xc*2.18947e-5)
2236
2237 xc=comexp(i,j,k,1)*comexp(i,j,k,1)*comexp(i,j,k,1)
2238 xc=xc*xc
2239 xc1=xc*xc
2240 xc=xc1*xc1
2241 xc=xc*xc
2242 comexp(i,j,k,2)=xc*xc1
2243
2244
2245
2246 else
2247
2248 xc=dcom(i,j,k)*(1.+(3.4268e-2+3.7401e-4*dt(i,j,k))*dt(i,j,k))
2249 comexp(i,j,k,1)=exp(-xc*4.74570e-5)
2250
2251 xc=comexp(i,j,k,1)*comexp(i,j,k,1)
2252 xc1=xc*xc
2253 xc2=xc1*xc1
2254 xc=xc2*xc2
2255 xc=xc*xc
2256 comexp(i,j,k,2)=xc1*xc2*xc
2257
2258 endif
2259
2260 enddo
2261 enddo
2262 enddo
2263
2264 return
2265 end
2266
2267 subroutine cfcexps(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,dcfc,dt,cfcexp)
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
2285
2286 implicit none
2287 integer ib,m,n,np,i,j,k
2288
2289
2290
a456aa407c Andr*2291 _RL dcfc(m,n,np),dt(m,n,np)
1662f365b2 Andr*2292
2293
2294
a456aa407c Andr*2295 _RL cfcexp(m,n,np)
1662f365b2 Andr*2296
2297
2298
a456aa407c Andr*2299 _RL a1,b1,fk1,a2,b2,fk2
1662f365b2 Andr*2300
2301
2302
a456aa407c Andr*2303 _RL xf
1662f365b2 Andr*2304
2305
2306
2307 do k=1,np
2308 do j=1,n
2309 do i=1,m
2310
2311
2312
2313 if (ib.eq.4) then
2314 xf=dcfc(i,j,k)*(1.+(a1+b1*dt(i,j,k))*dt(i,j,k))
2315 cfcexp(i,j,k)=exp(-xf*fk1)
2316 else
2317 xf=dcfc(i,j,k)*(1.+(a2+b2*dt(i,j,k))*dt(i,j,k))
2318 cfcexp(i,j,k)=exp(-xf*fk2)
2319 endif
2320
2321 enddo
2322 enddo
2323 enddo
2324
2325 return
2326 end
2327
2328 subroutine b10exps(m,n,np,dh2o,dcont,dco2,dn2o,pa,dt
2329 * ,h2oexp,conexp,co2exp,n2oexp)
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348 implicit none
2349 integer m,n,np,i,j,k
2350
2351
2352
a456aa407c Andr*2353 _RL dh2o(m,n,np),dcont(m,n,np),dn2o(m,n,np)
2354 _RL dco2(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*2355
2356
2357
a456aa407c Andr*2358 _RL h2oexp(m,n,np,6),conexp(m,n,np,3),co2exp(m,n,np,6,2)
1662f365b2 Andr*2359 * ,n2oexp(m,n,np,4)
2360
2361
2362
a456aa407c Andr*2363 _RL xx,xx1,xx2,xx3
1662f365b2 Andr*2364
2365
2366
2367 do k=1,np
2368 do j=1,n
2369 do i=1,m
2370
2371
2372
2373
2374 xx=dh2o(i,j,k)*(pa(i,j,k)/500.0)
2375 1 *(1.+(0.0149+6.20e-5*dt(i,j,k))*dt(i,j,k))
2376
2377
2378
2379
2380 h2oexp(i,j,k,1)=exp(-xx*0.10624)
2381
2382 xx=h2oexp(i,j,k,1)*h2oexp(i,j,k,1)
2383 xx=xx*xx
2384 h2oexp(i,j,k,2)=xx*xx
2385
2386 xx=h2oexp(i,j,k,2)*h2oexp(i,j,k,2)
2387 xx=xx*xx
2388 h2oexp(i,j,k,3)=xx*xx
2389
2390 xx=h2oexp(i,j,k,3)*h2oexp(i,j,k,3)
2391 xx=xx*xx
2392 h2oexp(i,j,k,4)=xx*xx
2393
2394 xx=h2oexp(i,j,k,4)*h2oexp(i,j,k,4)
2395 xx=xx*xx
2396 h2oexp(i,j,k,5)=xx*xx
2397
2398 xx=h2oexp(i,j,k,5)*h2oexp(i,j,k,5)
2399 xx=xx*xx
2400 h2oexp(i,j,k,6)=xx*xx
2401
2402
2403
2404
2405 xx=dco2(i,j,k)*(pa(i,j,k)/300.0)**0.5
2406 1 *(1.+(0.0179+1.02e-4*dt(i,j,k))*dt(i,j,k))
2407
2408
2409
2410
2411 co2exp(i,j,k,1,1)=exp(-xx*2.656e-5)
2412
2413 xx=co2exp(i,j,k,1,1)*co2exp(i,j,k,1,1)
2414 xx=xx*xx
2415 co2exp(i,j,k,2,1)=xx*xx
2416
2417 xx=co2exp(i,j,k,2,1)*co2exp(i,j,k,2,1)
2418 xx=xx*xx
2419 co2exp(i,j,k,3,1)=xx*xx
2420
2421 xx=co2exp(i,j,k,3,1)*co2exp(i,j,k,3,1)
2422 xx=xx*xx
2423 co2exp(i,j,k,4,1)=xx*xx
2424
2425 xx=co2exp(i,j,k,4,1)*co2exp(i,j,k,4,1)
2426 xx=xx*xx
2427 co2exp(i,j,k,5,1)=xx*xx
2428
2429 xx=co2exp(i,j,k,5,1)*co2exp(i,j,k,5,1)
2430 xx=xx*xx
2431 co2exp(i,j,k,6,1)=xx*xx
2432
2433
2434
2435
2436
2437 conexp(i,j,k,1)=exp(-dcont(i,j,k)*1.04995e+2)
2438
2439
2440
2441 xx=dn2o(i,j,k)*(1.+(1.4476e-3+3.6656e-6*dt(i,j,k))*dt(i,j,k))
2442
2443
2444
2445 n2oexp(i,j,k,1)=exp(-xx*0.25238)
2446
2447 xx=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)
2448 xx1=xx*xx
2449 xx1=xx1*xx1
2450 xx2=xx1*xx1
2451 xx3=xx2*xx2
2452 n2oexp(i,j,k,2)=xx*xx1*xx2*xx3
2453
2454 enddo
2455 enddo
2456 enddo
2457
2458 return
2459 end
2460
2461 subroutine tablup(k1,k2,m,n,np,nx,nh,nt,sabs,spre,stem,w1,p1,
2462 * dwe,dpe,coef1,coef2,coef3,tran)
2463
2464
2465
2466
2467
2468
2469
2470
2471
2472
2473
2474
2475
2476
2477
2478
2479
2480
2481
2482
2483
2484
2485
2486
2487
2488
2489
2490
2491
2492
2493
2494
2495
2496
2497
2498
2499 implicit none
2500 integer k1,k2,m,n,np,nx,nh,nt,i,j
2501
2502
2503
a456aa407c Andr*2504 _RL w1,p1,dwe,dpe
2505 _RL sabs(m,n,np+1),spre(m,n,np+1),stem(m,n,np+1)
2506 _RL coef1(nx,nh,nt),coef2(nx,nh,nt),coef3(nx,nh,nt)
1662f365b2 Andr*2507
2508
2509
a456aa407c Andr*2510 _RL tran(m,n)
1662f365b2 Andr*2511
2512
2513
a456aa407c Andr*2514 _RL x1,x2,x3,we,pe,fw,fp,pa,pb,pc,ax,ba,bb,t1,ca,cb,t2
1662f365b2 Andr*2515 integer iw,ip,nn
2516
2517
2518
2519 do j=1,n
2520 do i=1,m
2521
2522 nn=mod(i,nt)+1
2523
2524 x1=sabs(i,j,k2)-sabs(i,j,k1)
2525 x2=(spre(i,j,k2)-spre(i,j,k1))/x1
2526 x3=(stem(i,j,k2)-stem(i,j,k1))/x1
2527
2528 we=(log10(x1)-w1)/dwe
2529 pe=(log10(x2)-p1)/dpe
2530
2531 we=max(we,w1-2.*dwe)
2532 pe=max(pe,p1)
2533
2534 iw=int(we+1.5)
2535 ip=int(pe+1.5)
2536
2537 iw=min(iw,nh-1)
2538 iw=max(iw, 2)
2539
2540 ip=min(ip,nx-1)
2541 ip=max(ip, 1)
2542
2543 fw=we-float(iw-1)
2544 fp=pe-float(ip-1)
2545
2546
2547
2548 pa = coef1(ip,iw-1,nn)*(1.-fp)+coef1(ip+1,iw-1,nn)*fp
2549 pb = coef1(ip,iw, nn)*(1.-fp)+coef1(ip+1,iw, nn)*fp
2550 pc = coef1(ip,iw+1,nn)*(1.-fp)+coef1(ip+1,iw+1,nn)*fp
2551
2552
2553
2554 ax = (-pa*(1.-fw)+pc*(1.+fw)) *fw*0.5 + pb*(1.-fw*fw)
2555
2556
2557
2558 ba = coef2(ip,iw, nn)*(1.-fp)+coef2(ip+1,iw, nn)*fp
2559 bb = coef2(ip,iw+1,nn)*(1.-fp)+coef2(ip+1,iw+1,nn)*fp
2560 t1 = ba*(1.-fw) + bb*fw
2561
2562 ca = coef3(ip,iw, nn)*(1.-fp)+coef3(ip+1,iw, nn)*fp
2563 cb = coef3(ip,iw+1,nn)*(1.-fp)+coef3(ip+1,iw+1,nn)*fp
2564 t2 = ca*(1.-fw) + cb*fw
2565
2566
2567
2568 tran(i,j)= (ax + (t1+t2*x3) * x3)*tran(i,j)
2569
2570 enddo
2571 enddo
2572
2573 return
2574 end
2575
2576 subroutine h2okdis(ib,m,n,np,k,fkw,gkw,ne,h2oexp,conexp,
2577 * th2o,tcon,tran)
2578
2579
2580
2581
2582
2583
2584
2585
2586
2587
2588
2589
2590
2591
2592
2593
2594
2595
2596
2597
2598
2599
2600
2601
2602
2603
2604
2605
2606
2607 implicit none
2608 integer ib,m,n,np,k,i,j
2609
2610
2611
a456aa407c Andr*2612 _RL conexp(m,n,np,3),h2oexp(m,n,np,6)
1662f365b2 Andr*2613 integer ne(9)
a456aa407c Andr*2614 _RL fkw(6,9),gkw(6,3)
1662f365b2 Andr*2615
2616
2617
a456aa407c Andr*2618 _RL th2o(m,n,6),tcon(m,n,3),tran(m,n)
1662f365b2 Andr*2619
2620
2621
a456aa407c Andr*2622 _RL trnth2o
1662f365b2 Andr*2623
2624
2625
2626
2627
2628
2629
2630
2631
2632
2633
2634
2635
2636 do j=1,n
2637 do i=1,m
2638 th2o(i,j,1) = th2o(i,j,1)*h2oexp(i,j,k,1)
2639 th2o(i,j,2) = th2o(i,j,2)*h2oexp(i,j,k,2)
2640 th2o(i,j,3) = th2o(i,j,3)*h2oexp(i,j,k,3)
2641 th2o(i,j,4) = th2o(i,j,4)*h2oexp(i,j,k,4)
2642 th2o(i,j,5) = th2o(i,j,5)*h2oexp(i,j,k,5)
2643 th2o(i,j,6) = th2o(i,j,6)*h2oexp(i,j,k,6)
2644 enddo
2645 enddo
2646
2647 if (ne(ib).eq.0) then
2648
2649 do j=1,n
2650 do i=1,m
2651
2652 trnth2o =(fkw(1,ib)*th2o(i,j,1)
2653 * + fkw(2,ib)*th2o(i,j,2)
2654 * + fkw(3,ib)*th2o(i,j,3)
2655 * + fkw(4,ib)*th2o(i,j,4)
2656 * + fkw(5,ib)*th2o(i,j,5)
2657 * + fkw(6,ib)*th2o(i,j,6))
2658
2659 tran(i,j)=tran(i,j)*trnth2o
2660
2661 enddo
2662 enddo
2663
2664 elseif (ne(ib).eq.1) then
2665
2666 do j=1,n
2667 do i=1,m
2668
2669 tcon(i,j,1)= tcon(i,j,1)*conexp(i,j,k,1)
2670
2671 trnth2o =(fkw(1,ib)*th2o(i,j,1)
2672 * + fkw(2,ib)*th2o(i,j,2)
2673 * + fkw(3,ib)*th2o(i,j,3)
2674 * + fkw(4,ib)*th2o(i,j,4)
2675 * + fkw(5,ib)*th2o(i,j,5)
2676 * + fkw(6,ib)*th2o(i,j,6))*tcon(i,j,1)
2677
2678 tran(i,j)=tran(i,j)*trnth2o
2679
2680 enddo
2681 enddo
2682
2683 else
2684
2685 do j=1,n
2686 do i=1,m
2687
2688 tcon(i,j,1)= tcon(i,j,1)*conexp(i,j,k,1)
2689 tcon(i,j,2)= tcon(i,j,2)*conexp(i,j,k,2)
2690 tcon(i,j,3)= tcon(i,j,3)*conexp(i,j,k,3)
2691
2692
2693 trnth2o = ( gkw(1,1)*th2o(i,j,1)
2694 * + gkw(2,1)*th2o(i,j,2)
2695 * + gkw(3,1)*th2o(i,j,3)
2696 * + gkw(4,1)*th2o(i,j,4)
2697 * + gkw(5,1)*th2o(i,j,5)
2698 * + gkw(6,1)*th2o(i,j,6) ) * tcon(i,j,1)
2699 * + ( gkw(1,2)*th2o(i,j,1)
2700 * + gkw(2,2)*th2o(i,j,2)
2701 * + gkw(3,2)*th2o(i,j,3)
2702 * + gkw(4,2)*th2o(i,j,4)
2703 * + gkw(5,2)*th2o(i,j,5)
2704 * + gkw(6,2)*th2o(i,j,6) ) * tcon(i,j,2)
2705 * + ( gkw(1,3)*th2o(i,j,1)
2706 * + gkw(2,3)*th2o(i,j,2)
2707 * + gkw(3,3)*th2o(i,j,3)
2708 * + gkw(4,3)*th2o(i,j,4)
2709 * + gkw(5,3)*th2o(i,j,5)
2710 * + gkw(6,3)*th2o(i,j,6) ) * tcon(i,j,3)
2711
2712 tran(i,j)=tran(i,j)*trnth2o
2713
2714 enddo
2715 enddo
2716
2717 endif
2718
2719 return
2720 end
2721
2722 subroutine co2kdis(m,n,np,k,co2exp,tco2,tran)
2723
2724
2725
2726
2727
2728
2729
2730
2731
2732
2733
2734
2735
2736
2737
2738
2739
2740
2741 implicit none
2742 integer m,n,np,k,i,j
2743
2744
2745
a456aa407c Andr*2746 _RL co2exp(m,n,np,6,2)
1662f365b2 Andr*2747
2748
2749
a456aa407c Andr*2750 _RL tco2(m,n,6,2),tran(m,n)
1662f365b2 Andr*2751
2752
2753
a456aa407c Andr*2754 _RL xc
1662f365b2 Andr*2755
2756
2757
2758
2759 do j=1,n
2760 do i=1,m
2761
2762
2763
2764 tco2(i,j,1,1)=tco2(i,j,1,1)*co2exp(i,j,k,1,1)
2765 xc= 0.1395 *tco2(i,j,1,1)
2766
2767 tco2(i,j,2,1)=tco2(i,j,2,1)*co2exp(i,j,k,2,1)
2768 xc=xc+0.1407 *tco2(i,j,2,1)
2769
2770 tco2(i,j,3,1)=tco2(i,j,3,1)*co2exp(i,j,k,3,1)
2771 xc=xc+0.1549 *tco2(i,j,3,1)
2772
2773 tco2(i,j,4,1)=tco2(i,j,4,1)*co2exp(i,j,k,4,1)
2774 xc=xc+0.1357 *tco2(i,j,4,1)
2775
2776 tco2(i,j,5,1)=tco2(i,j,5,1)*co2exp(i,j,k,5,1)
2777 xc=xc+0.0182 *tco2(i,j,5,1)
2778
2779 tco2(i,j,6,1)=tco2(i,j,6,1)*co2exp(i,j,k,6,1)
2780 xc=xc+0.0220 *tco2(i,j,6,1)
2781
2782
2783
2784 tco2(i,j,1,2)=tco2(i,j,1,2)*co2exp(i,j,k,1,2)
2785 xc=xc+0.0766 *tco2(i,j,1,2)
2786
2787 tco2(i,j,2,2)=tco2(i,j,2,2)*co2exp(i,j,k,2,2)
2788 xc=xc+0.1372 *tco2(i,j,2,2)
2789
2790 tco2(i,j,3,2)=tco2(i,j,3,2)*co2exp(i,j,k,3,2)
2791 xc=xc+0.1189 *tco2(i,j,3,2)
2792
2793 tco2(i,j,4,2)=tco2(i,j,4,2)*co2exp(i,j,k,4,2)
2794 xc=xc+0.0335 *tco2(i,j,4,2)
2795
2796 tco2(i,j,5,2)=tco2(i,j,5,2)*co2exp(i,j,k,5,2)
2797 xc=xc+0.0169 *tco2(i,j,5,2)
2798
2799 tco2(i,j,6,2)=tco2(i,j,6,2)*co2exp(i,j,k,6,2)
2800 xc=xc+0.0059 *tco2(i,j,6,2)
2801
2802 tran(i,j)=tran(i,j)*xc
2803
2804 enddo
2805 enddo
2806
2807 return
2808 end
2809
2810 subroutine n2okdis(ib,m,n,np,k,n2oexp,tn2o,tran)
2811
2812
2813
2814
2815
2816
2817
2818
2819
2820
2821
2822
2823
2824
2825
2826
2827
2828
2829
2830 implicit none
2831 integer ib,m,n,np,k,i,j
2832
2833
2834
a456aa407c Andr*2835 _RL n2oexp(m,n,np,4)
1662f365b2 Andr*2836
2837
2838
a456aa407c Andr*2839 _RL tn2o(m,n,4),tran(m,n)
1662f365b2 Andr*2840
2841
2842
a456aa407c Andr*2843 _RL xc
1662f365b2 Andr*2844
2845
2846
2847
2848 do j=1,n
2849 do i=1,m
2850
2851
2852
2853 if (ib.eq.6) then
2854
2855 tn2o(i,j,1)=tn2o(i,j,1)*n2oexp(i,j,k,1)
2856 xc= 0.940414*tn2o(i,j,1)
2857
2858 tn2o(i,j,2)=tn2o(i,j,2)*n2oexp(i,j,k,2)
2859 xc=xc+0.059586*tn2o(i,j,2)
2860
2861
2862
2863 else
2864
2865 tn2o(i,j,1)=tn2o(i,j,1)*n2oexp(i,j,k,1)
2866 xc= 0.561961*tn2o(i,j,1)
2867
2868 tn2o(i,j,2)=tn2o(i,j,2)*n2oexp(i,j,k,2)
2869 xc=xc+0.138707*tn2o(i,j,2)
2870
2871 tn2o(i,j,3)=tn2o(i,j,3)*n2oexp(i,j,k,3)
2872 xc=xc+0.240670*tn2o(i,j,3)
2873
2874 tn2o(i,j,4)=tn2o(i,j,4)*n2oexp(i,j,k,4)
2875 xc=xc+0.058662*tn2o(i,j,4)
2876
2877 endif
2878
2879 tran(i,j)=tran(i,j)*xc
2880
2881 enddo
2882 enddo
2883
2884 return
2885 end
2886
2887 subroutine ch4kdis(ib,m,n,np,k,ch4exp,tch4,tran)
2888
2889
2890
2891
2892
2893
2894
2895
2896
2897
2898
2899
2900
2901
2902
2903
2904
2905
2906
2907 implicit none
2908 integer ib,m,n,np,k,i,j
2909
2910
2911
a456aa407c Andr*2912 _RL ch4exp(m,n,np,4)
1662f365b2 Andr*2913
2914
2915
a456aa407c Andr*2916 _RL tch4(m,n,4),tran(m,n)
1662f365b2 Andr*2917
2918
2919
a456aa407c Andr*2920 _RL xc
1662f365b2 Andr*2921
2922
2923
2924
2925 do j=1,n
2926 do i=1,m
2927
2928
2929
2930 if (ib.eq.6) then
2931
2932 tch4(i,j,1)=tch4(i,j,1)*ch4exp(i,j,k,1)
2933 xc= tch4(i,j,1)
2934
2935
2936
2937 else
2938
2939 tch4(i,j,1)=tch4(i,j,1)*ch4exp(i,j,k,1)
2940 xc= 0.610650*tch4(i,j,1)
2941
2942 tch4(i,j,2)=tch4(i,j,2)*ch4exp(i,j,k,2)
2943 xc=xc+0.280212*tch4(i,j,2)
2944
2945 tch4(i,j,3)=tch4(i,j,3)*ch4exp(i,j,k,3)
2946 xc=xc+0.107349*tch4(i,j,3)
2947
2948 tch4(i,j,4)=tch4(i,j,4)*ch4exp(i,j,k,4)
2949 xc=xc+0.001789*tch4(i,j,4)
2950
2951 endif
2952
2953 tran(i,j)=tran(i,j)*xc
2954
2955 enddo
2956 enddo
2957
2958 return
2959 end
2960
2961 subroutine comkdis(ib,m,n,np,k,comexp,tcom,tran)
2962
2963
2964
2965
2966
2967
2968
2969
2970
2971
2972
2973
2974
2975
2976
2977
2978
2979
2980
2981 implicit none
2982 integer ib,m,n,np,k,i,j
2983
2984
2985
a456aa407c Andr*2986 _RL comexp(m,n,np,2)
1662f365b2 Andr*2987
2988
2989
a456aa407c Andr*2990 _RL tcom(m,n,2),tran(m,n)
1662f365b2 Andr*2991
2992
2993
a456aa407c Andr*2994 _RL xc
1662f365b2 Andr*2995
2996
2997
2998
2999 do j=1,n
3000 do i=1,m
3001
3002
3003
3004 if (ib.eq.4) then
3005
3006 tcom(i,j,1)=tcom(i,j,1)*comexp(i,j,k,1)
3007 xc= 0.972025*tcom(i,j,1)
3008 tcom(i,j,2)=tcom(i,j,2)*comexp(i,j,k,2)
3009 xc=xc+0.027975*tcom(i,j,2)
3010
3011
3012
3013 else
3014
3015 tcom(i,j,1)=tcom(i,j,1)*comexp(i,j,k,1)
3016 xc= 0.961324*tcom(i,j,1)
3017 tcom(i,j,2)=tcom(i,j,2)*comexp(i,j,k,2)
3018 xc=xc+0.038676*tcom(i,j,2)
3019
3020 endif
3021
3022 tran(i,j)=tran(i,j)*xc
3023
3024 enddo
3025 enddo
3026
3027 return
3028 end
3029
3030 subroutine cfckdis(m,n,np,k,cfcexp,tcfc,tran)
3031
3032
3033
3034
3035
3036
3037
3038
3039
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049 implicit none
3050 integer m,n,np,k,i,j
3051
3052
3053
a456aa407c Andr*3054 _RL cfcexp(m,n,np)
1662f365b2 Andr*3055
3056
3057
a456aa407c Andr*3058 _RL tcfc(m,n),tran(m,n)
1662f365b2 Andr*3059
3060
3061
3062 do j=1,n
3063 do i=1,m
3064
3065 tcfc(i,j)=tcfc(i,j)*cfcexp(i,j,k)
3066 tran(i,j)=tran(i,j)*tcfc(i,j)
3067
3068 enddo
3069 enddo
3070
3071 return
3072 end
3073
3074 subroutine b10kdis(m,n,np,k,h2oexp,conexp,co2exp,n2oexp
3075 * ,th2o,tcon,tco2,tn2o,tran)
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105 implicit none
3106 integer m,n,np,k,i,j
3107
3108
3109
a456aa407c Andr*3110 _RL h2oexp(m,n,np,6),conexp(m,n,np,3),co2exp(m,n,np,6,2)
1662f365b2 Andr*3111 * ,n2oexp(m,n,np,4)
3112
3113
3114
a456aa407c Andr*3115 _RL th2o(m,n,6),tcon(m,n,3),tco2(m,n,6,2),tn2o(m,n,4)
1662f365b2 Andr*3116 * ,tran(m,n)
3117
3118
3119
a456aa407c Andr*3120 _RL xx
1662f365b2 Andr*3121
3122
3123
3124 do j=1,n
3125 do i=1,m
3126 tran(i,j)=1.0
3127 enddo
3128 enddo
3129
3130
3131
3132 do j=1,n
3133 do i=1,m
3134
3135 th2o(i,j,1)=th2o(i,j,1)*h2oexp(i,j,k,1)
3136 xx= 0.3153*th2o(i,j,1)
3137
3138 th2o(i,j,2)=th2o(i,j,2)*h2oexp(i,j,k,2)
3139 xx=xx+0.4604*th2o(i,j,2)
3140
3141 th2o(i,j,3)=th2o(i,j,3)*h2oexp(i,j,k,3)
3142 xx=xx+0.1326*th2o(i,j,3)
3143
3144 th2o(i,j,4)=th2o(i,j,4)*h2oexp(i,j,k,4)
3145 xx=xx+0.0798*th2o(i,j,4)
3146
3147 th2o(i,j,5)=th2o(i,j,5)*h2oexp(i,j,k,5)
3148 xx=xx+0.0119*th2o(i,j,5)
3149
3150 tran(i,j)=tran(i,j)*xx
3151
3152 enddo
3153 enddo
3154
3155
3156
3157 do j=1,n
3158 do i=1,m
3159
3160 tcon(i,j,1)=tcon(i,j,1)*conexp(i,j,k,1)
3161 tran(i,j)=tran(i,j)*tcon(i,j,1)
3162
3163 enddo
3164 enddo
3165
3166
3167
3168 do j=1,n
3169 do i=1,m
3170
3171 tco2(i,j,1,1)=tco2(i,j,1,1)*co2exp(i,j,k,1,1)
3172 xx= 0.2673*tco2(i,j,1,1)
3173
3174 tco2(i,j,2,1)=tco2(i,j,2,1)*co2exp(i,j,k,2,1)
3175 xx=xx+ 0.2201*tco2(i,j,2,1)
3176
3177 tco2(i,j,3,1)=tco2(i,j,3,1)*co2exp(i,j,k,3,1)
3178 xx=xx+ 0.2106*tco2(i,j,3,1)
3179
3180 tco2(i,j,4,1)=tco2(i,j,4,1)*co2exp(i,j,k,4,1)
3181 xx=xx+ 0.2409*tco2(i,j,4,1)
3182
3183 tco2(i,j,5,1)=tco2(i,j,5,1)*co2exp(i,j,k,5,1)
3184 xx=xx+ 0.0196*tco2(i,j,5,1)
3185
3186 tco2(i,j,6,1)=tco2(i,j,6,1)*co2exp(i,j,k,6,1)
3187 xx=xx+ 0.0415*tco2(i,j,6,1)
3188
3189 tran(i,j)=tran(i,j)*xx
3190
3191 enddo
3192 enddo
3193
3194
3195
3196 do j=1,n
3197 do i=1,m
3198
3199 tn2o(i,j,1)=tn2o(i,j,1)*n2oexp(i,j,k,1)
3200 xx= 0.970831*tn2o(i,j,1)
3201
3202 tn2o(i,j,2)=tn2o(i,j,2)*n2oexp(i,j,k,2)
3203 xx=xx+0.029169*tn2o(i,j,2)
3204
3205 tran(i,j)=tran(i,j)*(xx-1.0)
3206
3207 enddo
3208 enddo
3209
3210 return
3211 end