File indexing completed on 2023-03-03 06:10:02 UTC
view on githubraw file Latest commit 06d4643e on 2023-01-18 18:18:37 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
13fa5cb63c Jean*0002
0003
0004
e4ce4355da Jean*0005
0006
13fa5cb63c Jean*0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
e4ce4355da Jean*0024 SUBROUTINE TURBIO (im,jm,nlay,istrip,nymd,nhms,bi,bj, qbeg
6ab41f93a1 Andr*0025 1 ,ndturb,nltop,ptop, pz, uz, vz, tz, qz, ntracers,ptracers
1662f365b2 Andr*0026 2 ,plz,plze,dpres,pkht,pkz,ctmt,xxmt,yymt,zetamt,xlmt,khmt,tke
2a3ae9099b Andr*0027 3 ,tgz,fracland,landtype
1662f365b2 Andr*0028 4 ,tcanopy,ecanopy,tdeep,swetshal,swetroot,swetdeep,snodep,capac
c88fa11306 Andr*0029 5 ,nchp,nchptot,nchplnd,chfr,chlt,chlon,igrd,ityp,alai,agrn,thkz
0030 6 ,tprof
1662f365b2 Andr*0031 8 ,duturb, dvturb, dtturb,dqturb,radlwg,st4,dst4,radswg,radswt
0032 9 ,fdifpar,fdirpar,rainlsp,rainconv,snowfall,tempref
0033 1 ,imstturblw,imstturbsw,qliqavelw,qliqavesw,fccavelw,fccavesw
5cede9ba3f Andr*0034 2 ,qqgrid,myid)
5317312052 Jean*0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
e4ce4355da Jean*0045
5317312052 Jean*0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
1662f365b2 Andr*0111 implicit none
0112
6ab41f93a1 Andr*0113 integer im,jm,nlay,istrip,nymd,nhms,bi,bj,ndturb,nltop
175684e43e Andr*0114 integer ntracers, ptracers
c88fa11306 Andr*0115 integer nchp,nchptot,nchplnd
e4ce4355da Jean*0116 logical qbeg
a456aa407c Andr*0117 _RL ptop
5317312052 Jean*0118 _RL pz(im*jm,1),uz(im*jm,1,nlay),vz(im*jm,1,nlay)
0119 _RL tz(im*jm,1,nlay),qz(im*jm,1,nlay,ntracers)
0120 _RL plz(im*jm,1,nlay),plze(im*jm,1,nlay+1),dpres(im*jm,1,nlay)
0121 _RL pkht(im*jm,1,nlay+1),pkz(im*jm,1,nlay)
a456aa407c Andr*0122 _RL ctmt(nchp),xxmt(nchp),yymt(nchp),zetamt(nchp)
0123 _RL xlmt(nchp,nlay),khmt(nchp,nlay),tke(nchp,nlay)
5317312052 Jean*0124 _RL tgz(im*jm,1),fracland(im*jm,1)
0125 integer landtype(im*jm,1)
a456aa407c Andr*0126 _RL tcanopy(nchp),tdeep(nchp),ecanopy(nchp),swetshal(nchp)
0127 _RL swetroot(nchp),swetdeep(nchp),snodep(nchp),capac(nchp)
0128 _RL chfr(nchp),chlt(nchp),chlon(nchp)
5cede9ba3f Andr*0129 integer igrd(nchp),ityp(nchp)
5317312052 Jean*0130 _RL alai(nchp),agrn(nchp),thkz(im*jm,1)
1662f365b2 Andr*0131 logical tprof
5317312052 Jean*0132 _RL duturb(im*jm,1,nlay),dvturb(im*jm,1,nlay)
0133 _RL dtturb(im*jm,1,nlay),dqturb(im*jm,1,nlay,ntracers)
0134 _RL st4(im*jm,1),dst4(im*jm,1)
0135 _RL radswg(im*jm,1),radswt(im*jm,1),radlwg(im*jm,1)
0136 _RL fdifpar(im*jm,1),fdirpar(im*jm,1)
0137 _RL rainlsp(im*jm,1),rainconv(im*jm,1),snowfall(im*jm,1)
0138 _RL tempref (im*jm,1)
5cede9ba3f Andr*0139 integer imstturblw, imstturbsw
5317312052 Jean*0140 _RL qliqavesw(im*jm,1,nlay),qliqavelw(im*jm,1,nlay)
0141 _RL fccavelw (im*jm,1,nlay),fccavesw (im*jm,1,nlay)
0142 _RL qqgrid (im*jm,1,nlay)
5cede9ba3f Andr*0143 integer myid
0144
0145
1662f365b2 Andr*0146
0147 integer numstrips
0148 integer ijall
a456aa407c Andr*0149 _RL fmu,hice,tref,pref,cti,ed
1662f365b2 Andr*0150
5317312052 Jean*0151 parameter ( fmu = 0.00000 )
1662f365b2 Andr*0152 parameter ( hice = 300. )
5317312052 Jean*0153 parameter ( tref = 258. )
1662f365b2 Andr*0154 parameter ( pref = 500. )
0155 parameter ( cti = 0.0052 )
0156 parameter ( ed = 0.0 )
0157
5317312052 Jean*0158 _RL qliqtmp(im*jm,1,nlay)
0159 _RL fcctmp(im*jm,1,nlay)
0160 _RL tmpdiag(im*jm,1), tmpFac
a456aa407c Andr*0161 _RL thtgz(im*jm)
dee57c8146 Andr*0162 logical diagnostics_is_on
0163 external diagnostics_is_on
1662f365b2 Andr*0164
0165 integer nland
a456aa407c Andr*0166 _RL alwcoeff(nchp),blwcoeff(nchp)
0167 _RL netsw(nchp)
0168 _RL cnvprec(nchp),lsprec(nchp)
0169 _RL snowprec(nchp)
0170 _RL pardiff(nchp),pardirct(nchp)
0171 _RL pmsc(nchp)
0172 _RL netlw(nchp)
0173 _RL sqscat(nchp), rsoil1(nchp)
0174 _RL rsoil2(nchp)
0175 _RL rdc(nchp),u2fac(nchp)
0176 _RL z2ch(nchp)
0177 _RL zoch(nchp),cdrc(nchp)
0178 _RL cdsc(nchp)
0179 _RL dqsdt(nchp)
0180 _RL tground(nchp),qground(nchp)
0181 _RL utility(nchp)
0182 _RL qice(nchp)
0183 _RL dqice(nchp)
0184
0185 _RL dumsc(nchp,nlay),dvmsc(nchp,nlay)
0186 _RL dtmsc(nchp,nlay),dqmsc(nchp,nlay,ntracers)
0187
0188 _RL shg(nchp),z0(nchp),icethk(nchp)
1662f365b2 Andr*0189 integer water(nchp)
5317312052 Jean*0190
a456aa407c Andr*0191 _RL lats(istrip),lons(istrip),cosz(istrip),icest(istrip)
0192 _RL rainls(istrip),raincon(istrip),newsnow(istrip)
0193 _RL pardf(istrip),pardr(istrip),swnet(istrip)
daaf4a8d7f Andr*0194 _RL hlwdwn(istrip),alwrad(istrip),blwrad(istrip)
0195 _RL tmpnlay(istrip)
a456aa407c Andr*0196 _RL laistrip(istrip),grnstrip(istrip),z2str(istrip),cd(istrip)
0197 _RL scatstr(istrip), rs1str(istrip), rs2str(istrip)
0198 _RL rdcstr(istrip),u2fstr(istrip),dqsdtstr(istrip)
0199 _RL eturb(istrip),dedqa(istrip),dedtc(istrip)
0200 _RL hsturb(istrip),dhsdqa(istrip),dhsdtc(istrip)
0201 _RL savetc(istrip),saveqa(istrip),lwstrip(istrip)
0202 _RL chfrstr(istrip),psurf(istrip),shgstr(istrip)
1662f365b2 Andr*0203 integer types(istrip),igrdstr(istrip)
a456aa407c Andr*0204 _RL evap(istrip),shflux(istrip),runoff(istrip),bomb(istrip)
0205 _RL eint(istrip),esoi(istrip),eveg(istrip),esno(istrip)
0206 _RL smelt(istrip),hlatn(istrip),hlwup(istrip),gdrain(istrip)
0207 _RL runsrf(istrip),fwsoil(istrip),evpot(istrip)
0208 _RL strdg1(istrip),strdg2(istrip),strdg3(istrip),strdg4(istrip)
0209 _RL strdg5(istrip),strdg6(istrip),strdg7(istrip),strdg8(istrip)
0210 _RL strdg9(istrip),tmpstrip(istrip),qicestr(istrip)
0211 _RL dqicestr(istrip)
0212
0213 _RL u(istrip,nlay+1), v(istrip,nlay+1), th(istrip,nlay+1)
0214 _RL sh(istrip,nlay+1), thv(istrip,nlay+1), pe(istrip,nlay+1)
0215 _RL tracers(istrip,nlay+1,ntracers)
0216 _RL dpstr(istrip,nlay),pke(istrip,nlay+1)
0217 _RL pk(istrip,nlay), qq(istrip,nlay), p(istrip,nlay)
0218 _RL sri(istrip,nlay), skh(istrip,nlay), skm(istrip,nlay)
0219 _RL stuflux(istrip,nlay), stvflux(istrip,nlay)
0220 _RL sttflux(istrip,nlay), stqflux(istrip,nlay)
0221 _RL frqtrb(istrip,nlay-1)
0222 _RL dshdthg(istrip,nlay),dthdthg(istrip,nlay)
0223 _RL dshdshg(istrip,nlay),dthdshg(istrip,nlay)
0224 _RL transth(istrip,nlay), transsh(istrip,nlay)
0225
0226 _RL tc(istrip),td(istrip),qa(istrip)
0227 _RL swet1(istrip),swet2(istrip),swet3(istrip)
0228 _RL capacity(istrip),snowdepth(istrip)
0229 _RL stz0(istrip)
0230 _RL stdiag(istrip)
0231 _RL tends(istrip),sustar(istrip), sz0(istrip),pbldpth(istrip)
0232 _RL sct(istrip), scu(istrip), swinds(istrip)
0233 _RL stu2m(istrip),stv2m(istrip),stt2m(istrip),stq2m(istrip)
0234 _RL stu10m(istrip),stv10m(istrip),stt10m(istrip),stq10m(istrip)
1662f365b2 Andr*0235 integer stwatr(istrip)
a456aa407c Andr*0236 _RL wspeed(istrip)
1662f365b2 Andr*0237
daaf4a8d7f Andr*0238 _RL ctsave(istrip),xxsave(istrip),yysave(istrip)
0239 _RL zetasave(istrip)
a456aa407c Andr*0240 _RL xlsave(istrip,nlay),khsave(istrip,nlay)
0241 _RL qliq(istrip,nlay),turbfcc(istrip,nlay)
0242 _RL qliqmsc(nchp,nlay),fccmsc(nchp,nlay)
1662f365b2 Andr*0243
a456aa407c Andr*0244 _RL pi,secday,sdayopi2,rgas,akap,cp,alhl
0245 _RL faceps,grav,caltoj,virtcon,getcon
0246 _RL heatw,undef,timstp,delttrb,dttrb,ra
0247 _RL edle,rmu,cltj10,atimstp,tice,const
45ad04df92 Jean*0248 integer istnp1,istnlay,itrtrb,i,L,nn,nt
1662f365b2 Andr*0249 integer nocean, nice
45ad04df92 Jean*0250 integer ndmoist,time_left, ndum0,ndum1
1662f365b2 Andr*0251 integer ntracedim
e4ce4355da Jean*0252 _RL dtfac,timstp2
1662f365b2 Andr*0253
1a7d6fa3d1 Andr*0254 integer n,nsecf,nmonf,ndayf
0255 nsecf(n) = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)
0256 nmonf(n) = mod(n,10000)/100
0257 ndayf(n) = mod(n,100)
0258
06d4643e1f Jean*0259 #ifdef FIZHI_CRAY
4e1c053948 Jean*0260 #ifdef FIZHI_F77_COMPIL
1662f365b2 Andr*0261
0262 #endif
0263 #endif
0264
5317312052 Jean*0265
082e38725b Andr*0266
1662f365b2 Andr*0267 pi = 4.*atan(1.)
0268 secday = getcon('SDAY')
0269 sdayopi2 = getcon('SDAY') / (pi*2.)
0270 rgas = getcon('RGAS')
0271 akap = getcon('KAPPA')
0272 cp = getcon('CP')
0273 alhl = getcon('LATENT HEAT COND')
0274 faceps = getcon('EPSFAC')
0275 grav = getcon('GRAVITY')
0276 caltoj = getcon('CALTOJ')
0277 virtcon = getcon('VIRTCON')
0278 heatw = getcon('HEATW')
0279 undef = getcon('UNDEF')
0280 ntracedim= max(ntracers-ptracers,1)
0281
45ad04df92 Jean*0282 call get_alarm ( 'moist',ndum0,ndum1,ndmoist,time_left )
1662f365b2 Andr*0283 timstp = nsecf(ndturb)
0284 timstp2 = nsecf(ndmoist)
32362b8fa8 Cons*0285 dtfac = min( 1.0 _d 0, timstp/timstp2 )
1662f365b2 Andr*0286
5317312052 Jean*0287
0288
1662f365b2 Andr*0289 delttrb = nsecf(ndturb)
0290
0291 ijall = im * jm
0292 istnp1 = istrip * (nlay+1)
0293 istnlay = istrip * nlay
0294 itrtrb = ( timstp / delttrb ) + 0.1
0295 dttrb = timstp / float(itrtrb)
0296 edle = ed * 0.2
5317312052 Jean*0297
0298
13fa5cb63c Jean*0299
1662f365b2 Andr*0300 rmu = fmu * tref * rgas / pref
0301 cltj10 = 10. * caltoj
0302 atimstp = 1. / timstp
0303 tice = getcon('FREEZING-POINT')
0304
5317312052 Jean*0305
0306
0307
1662f365b2 Andr*0308
5317312052 Jean*0309
0310
0311 tmpFac = -1.
0312 CALL DIAGNOSTICS_SCALE_FILL( tgz,tmpFac,1,
0313 & 'DTG ',0,1,-3,bi,bj,myid)
0314
0315
0316
0317
0318
0319
1662f365b2 Andr*0320
c88fa11306 Andr*0321 numstrips = ((nchptot-1)/istrip) + 1
1662f365b2 Andr*0322
c88fa11306 Andr*0323 call grd2msc(pz(1,1),im,jm,igrd,pmsc,nchp,nchptot)
9524fe64b5 Andr*0324
c88fa11306 Andr*0325 call grd2msc(tgz,im,jm,igrd,tground,nchp,nchptot)
1662f365b2 Andr*0326 do i = 1,ijall
5317312052 Jean*0327 tmpdiag(i,1) = st4(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))
1662f365b2 Andr*0328 1 - dst4(i,1)* tgz(i,1)
0329 enddo
c88fa11306 Andr*0330 call grd2msc(tmpdiag,im,jm,igrd,alwcoeff,nchp,nchptot)
1662f365b2 Andr*0331 do i = 1,ijall
0332 tmpdiag(i,1) = dst4(i,1)
0333 enddo
c88fa11306 Andr*0334 call grd2msc(tmpdiag,im,jm,igrd,blwcoeff,nchp,nchptot)
1662f365b2 Andr*0335 do i = 1,ijall
0336 tmpdiag(i,1) = fdifpar(i,1) * radswt(i,1)
0337 enddo
c88fa11306 Andr*0338 call grd2msc(tmpdiag,im,jm,igrd,pardiff,nchp,nchptot)
1662f365b2 Andr*0339 do i = 1,ijall
0340 tmpdiag(i,1) = fdirpar(i,1) * radswt(i,1)
0341 enddo
c88fa11306 Andr*0342 call grd2msc(tmpdiag,im,jm,igrd,pardirct,nchp,nchptot)
1662f365b2 Andr*0343 do i = 1,ijall
0344 tmpdiag(i,1) = radswg(i,1) * radswt(i,1)
0345 enddo
c88fa11306 Andr*0346 call grd2msc(tmpdiag,im,jm,igrd,netsw,nchp,nchptot)
1662f365b2 Andr*0347 do i = 1,ijall
0348 tmpdiag(i,1) = radlwg(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))
0349 enddo
c88fa11306 Andr*0350 call grd2msc(tmpdiag,im,jm,igrd,netlw,nchp,nchptot)
0351 call grd2msc(thkz,im,jm,igrd,icethk,nchp,nchptot)
0352 call grd2msc(rainlsp,im,jm,igrd,lsprec,nchp,nchptot)
0353 call grd2msc(rainconv,im,jm,igrd,cnvprec,nchp,nchptot)
0354 call grd2msc(snowfall,im,jm,igrd,snowprec,nchp,nchptot)
1662f365b2 Andr*0355
0356
0357
0358 call chpprm(nymd,nhms,nchp,nchplnd,chlt,ityp,alai,
0359 1 agrn,zoch,z2ch,cdrc,cdsc,sqscat,u2fac,rsoil1,rsoil2,rdc)
0360
5317312052 Jean*0361
0362
0363
1662f365b2 Andr*0364
5317312052 Jean*0365
1662f365b2 Andr*0366
c88fa11306 Andr*0367 do i = 1,nchptot
1662f365b2 Andr*0368 water(i) = 0
0369 if((ityp(i).eq.100).and.(icethk(i).eq.0. ))water(i) = 1
0370 enddo
0371
5317312052 Jean*0372
13fa5cb63c Jean*0373
c88fa11306 Andr*0374 do i =1,nchptot
1662f365b2 Andr*0375 if (icethk(i).gt.0.) then
0376 z0(i) = 1.e-4
0377 else if (ityp(i).eq.100) then
0378 z0(i) = 3.e-4
0379 else
0380 z0(i) = zoch(i)
0381 endif
0382 enddo
0383
5317312052 Jean*0384
0385
1662f365b2 Andr*0386
0387 do i = 1,nchplnd
0388 tground(i) = tcanopy(i)
0389 enddo
0390
0391
0392
c88fa11306 Andr*0393 do I =1,nchptot
1662f365b2 Andr*0394 utility(I) = pmsc(i) + ptop
0395 call qsat ( tground(i),utility(i),shg(i),dqsdt(i),.true. )
0396 enddo
0397
5317312052 Jean*0398
0399
1662f365b2 Andr*0400
0401 do i = 1,nchplnd
0402 qground(i) = ecanopy(i)
0403 enddo
c88fa11306 Andr*0404 do i = nchplnd+1,nchptot
1662f365b2 Andr*0405 qground(i) = shg(i)
0406 enddo
0407
5317312052 Jean*0408
c88fa11306 Andr*0409 do i = nchplnd+1,nchptot
1662f365b2 Andr*0410 swetshal(i) = 1.
0411 enddo
0412
5317312052 Jean*0413
0414
1662f365b2 Andr*0415 const = ( cti / hice ) * cltj10
c88fa11306 Andr*0416 do i =1,nchptot
1662f365b2 Andr*0417 qice(i) = 0.0
0418 dqice(i) = 0.0
0419 if( icethk(i).gt.0.0 ) then
0420 qice(i) = const*(tice-tground(i))
0421 dqice(i) = -const
0422 endif
0423 enddo
0424
3a3c653dfd Andr*0425 if(diagnostics_is_on('QICE ',myid) ) then
0426 do i =1,ijall
0427 tmpdiag(i,1) = 0.0
0428 enddo
0429 call msc2grd (igrd,chfr,qice,nchp,nchptot,fracland,tmpdiag,im,jm)
0430 call diagnostics_fill(tmpdiag,'QICE ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0431 endif
0432
5317312052 Jean*0433
0434
0435
1662f365b2 Andr*0436
0437 do 2000 nn = 1, numstrips
0438
0439 call strip2tile(uz,igrd,u,nchp,ijall,istrip,nlay,nn)
0440 call strip2tile(vz,igrd,v,nchp,ijall,istrip,nlay,nn)
0441 call strip2tile(tz,igrd,th,nchp,ijall,istrip,nlay,nn)
0442 call strip2tile(qz(1,1,1,1),igrd,sh,nchp,ijall,istrip,nlay,nn)
0443 call strip2tile(dpres,igrd,dpstr,nchp,ijall,istrip,nlay,nn)
0444 call strip2tile(plz,igrd,p,nchp,ijall,istrip,nlay,nn)
0445 call strip2tile(plze,igrd,pe,nchp,ijall,istrip,nlay+1,nn)
0446 call strip2tile(pkz,igrd,pk,nchp,ijall,istrip,nlay,nn)
0447 call strip2tile(pkht,igrd,pke,nchp,ijall,istrip,nlay+1,nn)
257b288583 Andr*0448
0449
0450
0451
1662f365b2 Andr*0452
c88fa11306 Andr*0453 call stripit (z0,stz0,nchptot,nchp,istrip,1,nn)
0454 call stripit (tground,th(1,nlay+1),nchptot,nchp,istrip,1,nn)
0455 call stripit (pmsc,pe(1,nlay+1),nchptot,nchp,istrip,1,nn)
0456 call stripit (tke,qq,nchptot,nchp,istrip,nlay-1,nn)
0457 call stripit (ctmt,ctsave,nchptot,nchp,istrip,1,nn)
0458 call stripit (xxmt,xxsave,nchptot,nchp,istrip,1,nn)
0459 call stripit (yymt,yysave,nchptot,nchp,istrip,1,nn)
0460 call stripit (zetamt,zetasave,nchptot,nchp,istrip,1,nn)
0461 call stripit (xlmt,xlsave,nchptot,nchp,istrip,nlay,nn)
0462 call stripit (khmt,khsave,nchptot,nchp,istrip,nlay,nn)
0463 call stripitint (water,stwatr,nchptot,nchp,istrip,1,nn)
0464
0465 call stripitint (igrd,igrdstr,nchptot,nchp,istrip,1,nn)
0466 call stripit (chfr,chfrstr,nchptot,nchp,istrip,1,nn)
0467 call stripit (icethk,icest,nchptot,nchp,istrip,1,nn)
0468 call stripit (pardiff,pardf,nchptot,nchp,istrip,1,nn)
0469 call stripit (pardirct,pardr,nchptot,nchp,istrip,1,nn)
0470 call stripit (chlt,lats,nchptot,nchp,istrip,1,nn)
0471 call stripit (chlon,lons,nchptot,nchp,istrip,1,nn)
0472 call stripit (lsprec,rainls,nchptot,nchp,istrip,1,nn)
0473 call stripit (cnvprec,raincon,nchptot,nchp,istrip,1,nn)
0474 call stripit (snowprec,newsnow,nchptot,nchp,istrip,1,nn)
0475 call stripit (netsw,swnet,nchptot,nchp,istrip,1,nn)
0476 call stripit (netlw,lwstrip,nchptot,nchp,istrip,1,nn)
0477 call stripit (alwcoeff,alwrad,nchptot,nchp,istrip,1,nn)
0478 call stripit (blwcoeff,blwrad,nchptot,nchp,istrip,1,nn)
0479 call stripit (alai,laistrip,nchptot,nchp,istrip,1,nn)
0480 call stripit (agrn,grnstrip,nchptot,nchp,istrip,1,nn)
0481 call stripit (z2ch,z2str,nchptot,nchp,istrip,1,nn)
0482 call stripit (sqscat,scatstr,nchptot,nchp,istrip,1,nn)
0483 call stripit (rsoil1,rs1str,nchptot,nchp,istrip,1,nn)
0484 call stripit (rsoil2,rs2str,nchptot,nchp,istrip,1,nn)
0485 call stripit (rdc,rdcstr,nchptot,nchp,istrip,1,nn)
0486 call stripit (u2fac,u2fstr,nchptot,nchp,istrip,1,nn)
0487 call stripit (shg,shgstr,nchptot,nchp,istrip,1,nn)
0488 call stripit (dqsdt,dqsdtstr,nchptot,nchp,istrip,1,nn)
0489 call stripit ( qice, qicestr,nchptot,nchp,istrip,1,nn)
0490 call stripit (dqice,dqicestr,nchptot,nchp,istrip,1,nn)
0491 call stripitint (ityp,types,nchptot,nchp,istrip,1,nn)
0492
0493 call stripit (tground,tc,nchptot,nchp,istrip,1,nn)
0494 call stripit (tdeep,td,nchptot,nchp,istrip,1,nn)
0495 call stripit (qground,qa,nchptot,nchp,istrip,1,nn)
0496 call stripit (swetshal,swet1,nchptot,nchp,istrip,1,nn)
0497 call stripit (swetroot,swet2,nchptot,nchp,istrip,1,nn)
0498 call stripit (swetdeep,swet3,nchptot,nchp,istrip,1,nn)
0499 call stripit (snodep,snowdepth,nchptot,nchp,istrip,1,nn)
0500 call stripit (capac,capacity,nchptot,nchp,istrip,1,nn)
1662f365b2 Andr*0501
5317312052 Jean*0502 #ifdef FIZHI_USE_FIXED_DAY
0503 call astro ( 20040321,nhms,lats,lons,istrip,cosz,ra )
0504 #else
1662f365b2 Andr*0505 call astro ( nymd,nhms,lats,lons,istrip,cosz,ra )
5317312052 Jean*0506 #endif
1662f365b2 Andr*0507
5317312052 Jean*0508
1662f365b2 Andr*0509 nocean = 0
0510 nland = 0
0511 nice = 0
0512 do i = 1,istrip
0513 if( types(i).lt.100 ) nland = nland + 1
0514 if( types(i).eq.100 ) nocean = nocean + 1
0515 if( types(i).eq.100 .and. icest(i).gt.0.0 ) nice = nice + 1
0516 enddo
0517
5317312052 Jean*0518
0519
594697cf56 Andr*0520 do i =1,istrip
0521 u(i,nlay+1) = 0.
0522 v(i,nlay+1) = 0.
0523 enddo
0524
5317312052 Jean*0525
0526
1662f365b2 Andr*0527 do i =1,istrip
0528 th(i,nlay+1) = th(i,nlay+1) / pke(i,nlay+1)
0529 sh(i,nlay+1) = qa(i)
0530 enddo
0531
3a3c653dfd Andr*0532 if(diagnostics_is_on('QG ',myid) ) then
0533 do i=1,istrip
0534 tmpstrip(i) = sh(i,nlay+1)*1000
0535 enddo
0536 call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0537 & .false., 'QG ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0538 endif
1662f365b2 Andr*0539
5317312052 Jean*0540
0541
257b288583 Andr*0542
5317312052 Jean*0543
0544
0545
0546
1662f365b2 Andr*0547
5317312052 Jean*0548
0549
1662f365b2 Andr*0550 do L = 1,nlay+1
0551 do i =1,istrip
0552 thv(i,L) = 1. + virtcon * sh(i,L)
0553 thv(i,L) = th(i,L) * thv(i,L)
0554 enddo
0555 enddo
0556 do i =1,istrip
0557 sh(i,nlay+1) = qa(i)
0558 enddo
0559
5317312052 Jean*0560
1662f365b2 Andr*0561 do L =1,nlay
0562 do i =1,istrip
0563 qliq(i,L) = 0.
0564 turbfcc(i,L) = 0.
0565 enddo
0566 enddo
0567
5317312052 Jean*0568
0569
1662f365b2 Andr*0570 do i = 1,istrip
0571 eturb(i) = 0.
0572 scu(i) = 0.
0573 dedqa(i) = 0.
0574 dedtc(i) = 0.
0575 hsturb(i) = 0.
0576 dhsdqa(i) = 0.
0577 dhsdtc(i) = 0.
0578 enddo
0579
5317312052 Jean*0580
0581
3a3c653dfd Andr*0582 if(diagnostics_is_on('DTSRF ',myid) ) then
0583 do i=1,istrip
0584 stdiag(i) = ( thv(i,nlay+1)-thv(i,nlay) ) / pke(i,nlay+1)
0585 enddo
0586 call diag_vegtile_fill (stdiag,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0587 & .false., 'DTSRF ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0588 endif
257b288583 Andr*0589
5317312052 Jean*0590
0591
1662f365b2 Andr*0592 call trbflx(nn,th,thv,sh,u,v,qq,p,pe,pk,pke,dpstr,stwatr,stz0,
0593 1 tracers,ntracers-ptracers,ntracedim,dttrb,itrtrb,rmu,edle,qbeg,
0594 2 tprof,stuflux,stvflux,sri,skh,skm,swinds,sustar,sz0,frqtrb,
0595 3 pbldpth,sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,
d27a65ec88 Andr*0596 4 stq10m,istrip,nlay,nltop,nymd,nhms,grav,cp,rgas,faceps,virtcon,
0597 5 undef,dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,
1662f365b2 Andr*0598 6 hsturb,dhsdqa,dhsdtc,transth,transsh,
453e05dab3 Andr*0599 7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
1662f365b2 Andr*0600
c88fa11306 Andr*0601 call pastit (qq,tke,istrip,nchp,nchptot,nlay,nn)
0602 call pastit (ctsave,ctmt,istrip,nchp,nchptot,1,nn)
0603 call pastit (xxsave,xxmt,istrip,nchp,nchptot,1,nn)
0604 call pastit (yysave,yymt,istrip,nchp,nchptot,1,nn)
0605 call pastit (zetasave,zetamt,istrip,nchp,nchptot,1,nn)
0606 call pastit (xlsave,xlmt,istrip,nchp,nchptot,nlay,nn)
0607 call pastit (khsave,khmt,istrip,nchp,nchptot,nlay,nn)
1662f365b2 Andr*0608
c88fa11306 Andr*0609 call pastit (qliq ,qliqmsc,istrip,nchp,nchptot,nlay,nn)
0610 call pastit (turbfcc,fccmsc,istrip,nchp,nchptot,nlay,nn)
1662f365b2 Andr*0611
5317312052 Jean*0612
1662f365b2 Andr*0613 do i = 1,istrip
0614 evpot(i) = transsh(i,nlay) * (shgstr(i) - sh(i,nlay))
0615 enddo
0616
0617
0618
0619
5317312052 Jean*0620
1662f365b2 Andr*0621 do i = 1,istrip
0622 savetc(i) = tc(i)
0623 saveqa(i) = qa(i)
0624 enddo
0625 do i = 1,istrip
32362b8fa8 Cons*0626 cosz(i) = max(cosz(i),0.0001 _d 0)
1662f365b2 Andr*0627 cd(i) = scu(i)*scu(i)
0628 tmpnlay(i) = th(i,nlay)*pk(i,nlay)
0629 hlwdwn(i) = alwrad(i)+blwrad(i)*tc(i)-lwstrip(i)
0630 psurf(i) = pe(i,nlay+1)
0631 wspeed(i) = sqrt(u(i,nlay)*u(i,nlay) + v(i,nlay)*v(i,nlay))
44d19a8651 Andr*0632 if(wspeed(i) .lt. 1.e-20) wspeed(i) = 1.e-20
1662f365b2 Andr*0633
13fa5cb63c Jean*0634
0635
0636
1662f365b2 Andr*0637 enddo
0638
0639 do i = 1,istrip
0640 eturb(i) = eturb(i) * pke(i,nlay+1)
0641 dedqa(i) = dedqa(i) * pke(i,nlay+1)
0642 hsturb(i) = hsturb(i) * pke(i,nlay+1)
0643 enddo
0644
0645 do i = 1,istrip
0646 strdg1(i) = 0.
0647 strdg2(i) = 0.
0648 strdg3(i) = 0.
0649 strdg4(i) = 0.
0650 strdg5(i) = 0.
0651 strdg6(i) = 0.
0652 strdg7(i) = 0.
0653 strdg8(i) = 0.
0654 strdg9(i) = 0.
0655 bomb(i) = 0.
0656 runoff(i) = 0.
0657 eint(i) = 0.
0658 esoi(i) = 0.
0659 eveg(i) = 0.
0660 esno(i) = 0.
0661 smelt(i) = 0.
0662 hlatn(i) = 0.
0663 hlwup(i) = 0.
0664 gdrain(i) = 0.
0665 runsrf(i) = 0.
0666 fwsoil(i) = 0.
0667 enddo
0668
5317312052 Jean*0669
0670
0671
3a3c653dfd Andr*0672 if(diagnostics_is_on('SNOWFALL',myid) ) then
0673 do i = 1,istrip
5317312052 Jean*0674 tmpstrip(i) = newsnow(i)*86400
3a3c653dfd Andr*0675 enddo
0676 call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0677 & .false., 'SNOWFALL', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0678 endif
dee57c8146 Andr*0679 if(diagnostics_is_on('RAINCON ',myid) ) then
0680 do i = 1,istrip
5317312052 Jean*0681 tmpstrip(i) = raincon(i)*86400
dee57c8146 Andr*0682 enddo
0683 call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0684 & .false., 'RAINCON ', 1, 1, bi, bj, myid)
dee57c8146 Andr*0685 endif
0686 if(diagnostics_is_on('RAINLSP ',myid) ) then
0687 do i = 1,istrip
5317312052 Jean*0688 tmpstrip(i) = rainls(i)*86400
dee57c8146 Andr*0689 enddo
0690 call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0691 & .false., 'RAINLSP ', 1, 1, bi, bj, myid)
dee57c8146 Andr*0692 endif
3a3c653dfd Andr*0693 if(diagnostics_is_on('GREEN ',myid) ) then
0694 call diag_vegtile_fill (grnstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0695 & .false., 'GREEN ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0696 endif
0697 if(diagnostics_is_on('LAI ',myid) ) then
0698 call diag_vegtile_fill (laistrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0699 & .false., 'LAI ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0700 endif
dee57c8146 Andr*0701 if(diagnostics_is_on('PARDR ',myid) ) then
0702 call diag_vegtile_fill (pardr,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0703 & .false., 'PARDR ', 1, 1, bi, bj, myid)
dee57c8146 Andr*0704 endif
3a3c653dfd Andr*0705 if(diagnostics_is_on('PARDF ',myid) ) then
0706 call diag_vegtile_fill (pardf,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0707 & .false., 'PARDF ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0708 endif
0709 if(diagnostics_is_on('DLWDTC ',myid) ) then
0710 call diag_vegtile_fill (blwrad,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0711 & .false., 'DLWDTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0712 endif
0713 if(diagnostics_is_on('DHDTC ',myid) ) then
0714 call diag_vegtile_fill (dhsdtc,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0715 & .false., 'DHDTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0716 endif
0717 if(diagnostics_is_on('DEDTC ',myid) ) then
0718 call diag_vegtile_fill (dedtc,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0719 & .false., 'DEDTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0720 endif
0721 if(diagnostics_is_on('DHDQA ',myid) ) then
0722 call diag_vegtile_fill (dhsdqa,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0723 & .false., 'DHDQA ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0724 endif
0725 if(diagnostics_is_on('DEDQA ',myid) ) then
0726 call diag_vegtile_fill (dedqa,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0727 & .false., 'DEDQA ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0728 endif
0729 if(diagnostics_is_on('LWGDOWN ',myid) ) then
0730 call diag_vegtile_fill (hlwdwn,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0731 & .false., 'LWGDOWN ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0732 endif
5317312052 Jean*0733
1662f365b2 Andr*0734
44d19a8651 Andr*0735 if(nland.gt.0)then
257b288583 Andr*0736
44d19a8651 Andr*0737 call tile (
0738 I nland, timstp, types, rainls, raincon, newsnow, wspeed,
0739 I eturb, dedqa, dedtc, hsturb, dhsdqa, dhsdtc,
0740 I tmpnlay, sh(1,nlay), cd, cosz, pardr, pardf,
0741 I swnet, hlwdwn, psurf, laistrip, grnstrip, z2str,
0742 I scatstr, rs1str, rs2str, rdcstr, u2fstr,
0743 I shgstr, dqsdtstr, alwrad, blwrad,
0744 U tc, td, qa, swet1, swet2, swet3, capacity, snowdepth,
0745 O evap, shflux, runoff, bomb,
0746 O eint, esoi, eveg, esno, smelt, hlatn,
0747 O hlwup, gdrain, runsrf, fwsoil,
0748 O strdg1, strdg2, strdg3, strdg4,
0749 O strdg5, strdg6, strdg7, strdg8, strdg9)
0750 endif
0751
0752 if( nice.gt.0 ) then
0753 call seaice ( nocean, timstp, hice,
13fa5cb63c Jean*0754 & eturb(nland+1), dedtc(nland+1),
0755 & hsturb(nland+1), dhsdtc(nland+1),
0756 & qicestr(nland+1), dqicestr(nland+1),
0757 & swnet(nland+1), lwstrip(nland+1), blwrad(nland+1),
0758 & pke(nland+1,nlay+1), icest(nland+1),
0759 & tc(nland+1), qa(nland+1) )
44d19a8651 Andr*0760 endif
1662f365b2 Andr*0761
5317312052 Jean*0762
0763
0764
1662f365b2 Andr*0765
3a3c653dfd Andr*0766 if(diagnostics_is_on('RUNOFF ',myid) ) then
0767 call diag_vegtile_fill (runoff,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0768 & .false., 'RUNOFF ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0769 endif
0770 if(diagnostics_is_on('FWSOIL ',myid) ) then
0771 call diag_vegtile_fill (fwsoil,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0772 & .false., 'FWSOIL ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0773 endif
0774 if(diagnostics_is_on('GDRAIN ',myid) ) then
0775 call diag_vegtile_fill (gdrain,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0776 & .false., 'GDRAIN ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0777 endif
0778 if(diagnostics_is_on('SNOWMELT',myid) ) then
0779 call diag_vegtile_fill (smelt,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0780 & .false., 'SNOWMELT', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0781 endif
0782 if(diagnostics_is_on('EVEG ',myid) ) then
0783 call diag_vegtile_fill (eveg,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0784 & .false., 'EVEG ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0785 endif
0786 if(diagnostics_is_on('ESNOW ',myid) ) then
c9f33fa462 Andr*0787 call diag_vegtile_fill (esno,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0788 & .false., 'ESNOW ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0789 endif
0790 if(diagnostics_is_on('ESOIL ',myid) ) then
c9f33fa462 Andr*0791 call diag_vegtile_fill (esoi,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0792 & .false., 'ESOIL ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0793 endif
0794 if(diagnostics_is_on('ERESV ',myid) ) then
0795 call diag_vegtile_fill (eint,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0796 & .false., 'ERESV ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0797 endif
0798 if(diagnostics_is_on('EVPOT ',myid) ) then
0799 call diag_vegtile_fill (evpot,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0800 & .false., 'EVPOT ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0801 endif
0802 if(diagnostics_is_on('DTC ',myid) ) then
0803 call diag_vegtile_fill (strdg1,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0804 & .false., 'DTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0805 endif
0806 if(diagnostics_is_on('DQC ',myid) ) then
0807 call diag_vegtile_fill (strdg2,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0808 & .false., 'DQC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0809 endif
0810 if(diagnostics_is_on('TCDTC ',myid) ) then
0811 call diag_vegtile_fill (strdg3,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0812 & .false., 'TCDTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0813 endif
0814 if(diagnostics_is_on('RADDTC ',myid) ) then
0815 call diag_vegtile_fill (strdg4,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0816 & .false., 'RADDTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0817 endif
0818 if(diagnostics_is_on('SENSDTC ',myid) ) then
0819 call diag_vegtile_fill (strdg5,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0820 & .false., 'SENSDTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0821 endif
0822 if(diagnostics_is_on('LATDTC ',myid) ) then
0823 call diag_vegtile_fill (strdg6,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0824 & .false., 'LATDTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0825 endif
0826 if(diagnostics_is_on('TDDTC ',myid) ) then
0827 call diag_vegtile_fill (strdg7,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0828 & .false., 'TDDTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0829 endif
0830 if(diagnostics_is_on('QCDTC ',myid) ) then
0831 call diag_vegtile_fill (strdg8,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0832 & .false., 'QCDTC ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0833 endif
5317312052 Jean*0834
1662f365b2 Andr*0835
c88fa11306 Andr*0836 call pastit (tc,tground,istrip,nchp,nchptot,1,nn)
0837 call pastit (td,tdeep,istrip,nchp,nchptot,1,nn)
0838 call pastit (qa,qground,istrip,nchp,nchptot,1,nn)
0839 call pastit (swet1,swetshal,istrip,nchp,nchptot,1,nn)
0840 call pastit (swet2,swetroot,istrip,nchp,nchptot,1,nn)
0841 call pastit (swet3,swetdeep,istrip,nchp,nchptot,1,nn)
0842 call pastit (capacity,capac,istrip,nchp,nchptot,1,nn)
0843 call pastit (snowdepth,snodep,istrip,nchp,nchptot,1,nn)
1662f365b2 Andr*0844
5317312052 Jean*0845
0846
0847
1662f365b2 Andr*0848
0849 do i =1,istrip
0850 th(i,nlay+1) = tc(i) / pke(i,nlay+1)
0851 enddo
0852 do L = 1,nlay
0853 do i =1,istrip
0854 th(i,L) = th(i,L) + dthdthg(i,L)*(tc(i)-savetc(i))/pke(i,nlay+1)
0855 enddo
0856 enddo
0857
0858 do i =1,istrip
0859 sh(i,nlay+1) = qa(i)
0860 enddo
0861 do L = 1,nlay
0862 do i =1,istrip
0863 sh(i,L) = sh(i,L) + dshdshg(i,L)*(qa(i)-saveqa(i))
0864 enddo
0865 enddo
0866
0867 do L = 1,nlay
0868 do i =1,istrip
0869 sttflux(i,L) = transth(i,L) * (th(i,L+1)-th(i,L))
0870 stqflux(i,L) = transsh(i,L) * (sh(i,L+1)-sh(i,L))
0871 enddo
0872 enddo
5317312052 Jean*0873
0874
0875
1662f365b2 Andr*0876 do l=1,nlay
0877 call strip2tile(uz(1,1,l),igrd,tmpstrip,nchp,ijall,
0878 1 istrip,1,nn)
0879 do i =1,istrip
0880 tends(i) = ( u(i,l)-tmpstrip(i) )
0881 enddo
c88fa11306 Andr*0882 call pastit (tends,dumsc(1,l),istrip,nchp,nchptot,1,nn)
1662f365b2 Andr*0883
0884 call strip2tile(vz(1,1,l),igrd,tmpstrip,nchp,ijall,
0885 1 istrip,1,nn)
0886 do i =1,istrip
0887 tends(i) = ( v(i,l)-tmpstrip(i) )
0888 enddo
c88fa11306 Andr*0889 call pastit (tends,dvmsc(1,l),istrip,nchp,nchptot,1,nn)
1662f365b2 Andr*0890
0891 call strip2tile(tz(1,1,l),igrd,tmpstrip,nchp,ijall,
0892 1 istrip,1,nn)
0893 do i =1,istrip
0894 tends(i) = ( th(i,l)-tmpstrip(i) )
0895 enddo
44d19a8651 Andr*0896
c88fa11306 Andr*0897 call pastit (tends,dtmsc(1,l),istrip,nchp,nchptot,1,nn)
1662f365b2 Andr*0898
0899 call strip2tile(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,
0900 1 istrip,1,nn)
0901 do i =1,istrip
0902 tends(i) = ( sh(i,l)-tmpstrip(i) )
0903 enddo
44d19a8651 Andr*0904
c88fa11306 Andr*0905 call pastit (tends,dqmsc(1,l,1),istrip,nchp,nchptot,1,nn)
1662f365b2 Andr*0906
257b288583 Andr*0907
0908
0909
0910
0911
0912
0913
0914
0915
1662f365b2 Andr*0916
0917 enddo
5317312052 Jean*0918
0919
0920
1662f365b2 Andr*0921
5317312052 Jean*0922
0923
0924
1662f365b2 Andr*0925
3a3c653dfd Andr*0926 if(diagnostics_is_on('EVAP ',myid) ) then
0927 do i=1,istrip
0928 tmpstrip(i) = stqflux(i,nlay) * 86400
0929 enddo
0930 call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0931 & .false., 'EVAP ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0932 endif
0933 if(diagnostics_is_on('EFLUX ',myid) ) then
0934 do i=1,istrip
0935 tmpstrip(i) = stqflux(i,nlay) * alhl
0936 enddo
0937 call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0938 & .false., 'EFLUX ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0939 endif
0940 if(diagnostics_is_on('HFLUX ',myid) ) then
0941 do i=1,istrip
0942 tmpstrip(i) = sttflux(i,nlay) * cp
0943 enddo
0944 call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0945 & .false., 'HFLUX ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0946 endif
a01198c294 Andr*0947 if(diagnostics_is_on('TUFLUX ',myid) ) then
0948 call diag_vegtile_fill (stuflux,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0949 & .false., 'TUFLUX ', 0, nlay, bi, bj, myid)
a01198c294 Andr*0950 endif
0951 if(diagnostics_is_on('TVFLUX ',myid) ) then
0952 call diag_vegtile_fill (stvflux,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0953 & .false., 'TVFLUX ', 0, nlay, bi, bj, myid)
a01198c294 Andr*0954 endif
0955 if(diagnostics_is_on('TTFLUX ',myid) ) then
0956 do l=1,nlay
0957 do i=1,istrip
0958 sttflux(i,l) = sttflux(i,l) * cp
0959 enddo
0960 enddo
0961 call diag_vegtile_fill (sttflux,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0962 & .false., 'TTFLUX ', 0, nlay, bi, bj, myid)
a01198c294 Andr*0963 endif
0964 if(diagnostics_is_on('TQFLUX ',myid) ) then
0965 do l=1,nlay
0966 do i=1,istrip
0967 stqflux(i,l) = stqflux(i,l) * alhl
0968 enddo
0969 enddo
0970 call diag_vegtile_fill (stqflux,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0971 & .false., 'TQFLUX ', 0, nlay, bi, bj, myid)
a01198c294 Andr*0972 endif
dee57c8146 Andr*0973 if(diagnostics_is_on('RI ',myid) ) then
0974 call diag_vegtile_fill (sri,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0975 & .false., 'RI ', 0, nlay, bi, bj, myid)
dee57c8146 Andr*0976 endif
0977 if(diagnostics_is_on('KH ',myid) ) then
0978 call diag_vegtile_fill (skh,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0979 & .false., 'KH ', 0, nlay, bi, bj, myid)
dee57c8146 Andr*0980 endif
0981 if(diagnostics_is_on('KM ',myid) ) then
0982 call diag_vegtile_fill (skm,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0983 & .false., 'KM ', 0, nlay, bi, bj, myid)
dee57c8146 Andr*0984 endif
3a3c653dfd Andr*0985 if(diagnostics_is_on('CT ',myid) ) then
0986 call diag_vegtile_fill (sct,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0987 & .false., 'CT ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0988 endif
0989 if(diagnostics_is_on('CU ',myid) ) then
0990 call diag_vegtile_fill (scu,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0991 & .false., 'CU ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0992 endif
0993 if(diagnostics_is_on('WINDS ',myid) ) then
0994 call diag_vegtile_fill (swinds,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*0995 & .false., 'WINDS ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*0996 endif
0997 if(diagnostics_is_on('UFLUX ',myid) ) then
0998 call diag_vegtile_fill (stuflux(1,nlay),igrd,chfrstr,istrip,nchp,
13fa5cb63c Jean*0999 & nn,.false., 'UFLUX ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1000 endif
1001 if(diagnostics_is_on('VFLUX ',myid) ) then
1002 call diag_vegtile_fill (stvflux(1,nlay),igrd,chfrstr,istrip,nchp,
13fa5cb63c Jean*1003 & nn,.false., 'VFLUX ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1004 endif
1005 if(diagnostics_is_on('USTAR ',myid) ) then
1006 call diag_vegtile_fill (sustar,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1007 & .false., 'USTAR ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1008 endif
1009 if(diagnostics_is_on('Z0 ',myid) ) then
1010 call diag_vegtile_fill (sz0,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1011 & .false., 'Z0 ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1012 endif
1013 if(diagnostics_is_on('FRQTRB ',myid) ) then
1014 call diag_vegtile_fill (frqtrb,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1015 & .false., 'FRQTRB ', 0, nlay-1, bi, bj, myid)
3a3c653dfd Andr*1016 endif
1017 if(diagnostics_is_on('PBL ',myid) ) then
1018 call diag_vegtile_fill (pbldpth,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1019 & .false., 'PBL ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1020 endif
1021 if(diagnostics_is_on('U2M ',myid) ) then
1022 call diag_vegtile_fill (stu2m,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1023 & .false., 'U2M ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1024 endif
1025 if(diagnostics_is_on('V2M ',myid) ) then
1026 call diag_vegtile_fill (stv2m,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1027 & .false., 'V2M ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1028 endif
1029 if(diagnostics_is_on('T2M ',myid) ) then
1030 call diag_vegtile_fill (stt2m,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1031 & .false., 'T2M ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1032 endif
1033 if(diagnostics_is_on('Q2M ',myid) ) then
1034 do i=1,istrip
1035 if( stq2m(i).ne.undef ) then
1036 tmpstrip(i) = stq2m(i) * 1000
1037 else
1038 tmpstrip(i) = undef
1039 endif
1040 enddo
1041 call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1042 & .false., 'Q2M ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1043 endif
1044 if(diagnostics_is_on('U10M ',myid) ) then
1045 call diag_vegtile_fill (stu10m,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1046 & .false., 'U10M ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1047 endif
1048 if(diagnostics_is_on('V10M ',myid) ) then
1049 call diag_vegtile_fill (stv10m,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1050 & .false., 'V10M ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1051 endif
1052 if(diagnostics_is_on('T10M ',myid) ) then
1053 call diag_vegtile_fill (stt10m,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1054 & .false., 'T10M ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1055 endif
1056 if(diagnostics_is_on('Q10M ',myid) ) then
1057 do i=1,istrip
1058 if( stq10m(i).ne.undef ) then
1059 tmpstrip(i) = stq10m(i) * 1000
1060 else
1061 tmpstrip(i) = undef
1062 endif
1063 enddo
1064 call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1065 & .false., 'Q10M ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1066 endif
1662f365b2 Andr*1067
5317312052 Jean*1068
1069
1070
1071
3a3c653dfd Andr*1072 if(diagnostics_is_on('TDEEP ',myid) ) then
1073 call diag_vegtile_fill (td,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1074 & .false., 'TDEEP ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1075 endif
1076 if(diagnostics_is_on('QCANOPY ',myid) ) then
1077 call diag_vegtile_fill (qa,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1078 & .false., 'QCANOPY ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1079 endif
1080 if(diagnostics_is_on('SMSHAL ',myid) ) then
1081 call diag_vegtile_fill (swet1,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1082 & .false., 'SMSHAL ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1083 endif
1084 if(diagnostics_is_on('SMROOT ',myid) ) then
1085 call diag_vegtile_fill (swet2,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1086 & .false., 'SMROOT ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1087 endif
1088 if(diagnostics_is_on('SMDEEP ',myid) ) then
1089 call diag_vegtile_fill (swet3,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1090 & .false., 'SMDEEP ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1091 endif
1092 if(diagnostics_is_on('CAPACITY',myid) ) then
1093 call diag_vegtile_fill (capacity,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1094 & .false., 'CAPACITY', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1095 endif
1096 if(diagnostics_is_on('SNOW ',myid) ) then
1097 call diag_vegtile_fill (snowdepth,igrd,chfrstr,istrip,nchp,nn,
13fa5cb63c Jean*1098 & .false., 'SNOW ', 1, 1, bi, bj, myid)
3a3c653dfd Andr*1099 endif
5317312052 Jean*1100
1101
1662f365b2 Andr*1102
1103 2000 continue
1104
5317312052 Jean*1105
1106
1107
1108
1662f365b2 Andr*1109 imstturblw = imstturblw + 1
1110 imstturbsw = imstturbsw + 1
1111
5317312052 Jean*1112
1113
c88fa11306 Andr*1114 do i =1,nchptot
1662f365b2 Andr*1115 if( (icethk(i).gt.0.).and.(tground(i).gt.tice) ) tground(i)=tice
1116 enddo
1117
5317312052 Jean*1118
1119
1120
3f83b87ea3 Andr*1121 do i =1,nchptot
1662f365b2 Andr*1122 tcanopy(i) = tground(i)
1123 ecanopy(i) = qground(i)
1124 enddo
1125
1126
5317312052 Jean*1127
1662f365b2 Andr*1128 do L = 1,nlay
1129 do i = 1,ijall
1130 duturb(i,1,L) = 0.
1131 dvturb(i,1,L) = 0.
1132 dtturb(i,1,L) = 0.
1133 qqgrid(i,1,L) = 0.
1134 qliqtmp(i,1,L) = 0.
1135 fcctmp(i,1,L) = 0.
1136 enddo
1137 do nt = 1,ntracers
1138 do i = 1,ijall
1139 dqturb(i,1,L,nt) = 0.
1140 enddo
1141 enddo
1142 enddo
1143
1144
5317312052 Jean*1145
1662f365b2 Andr*1146 do l = 1,nlay
c88fa11306 Andr*1147 call msc2grd(igrd,chfr,dumsc(1,L),nchp,nchptot,fracland,
13fa5cb63c Jean*1148 & duturb(1,1,L),im,jm)
c88fa11306 Andr*1149 call msc2grd(igrd,chfr,dvmsc(1,L),nchp,nchptot,fracland,
13fa5cb63c Jean*1150 & dvturb(1,1,L),im,jm)
c88fa11306 Andr*1151 call msc2grd(igrd,chfr,dtmsc(1,L),nchp,nchptot,fracland,
13fa5cb63c Jean*1152 & dtturb(1,1,L),im,jm)
1662f365b2 Andr*1153 do nt = 1,ntracers
c88fa11306 Andr*1154 call msc2grd(igrd,chfr,dqmsc(1,L,nt),nchp,nchptot,fracland,
13fa5cb63c Jean*1155 & dqturb(1,1,L,nt),im,jm)
1662f365b2 Andr*1156 enddo
c88fa11306 Andr*1157 call msc2grd(igrd,chfr, tke(1,L),nchp,nchptot,fracland,
13fa5cb63c Jean*1158 & qqgrid(1,1,L),im,jm)
1662f365b2 Andr*1159
5317312052 Jean*1160 call msc2grd(igrd,chfr, fccmsc(1,L),nchp,nchptot,fracland,
13fa5cb63c Jean*1161 & fcctmp(1,1,L),im,jm)
c88fa11306 Andr*1162 call msc2grd(igrd,chfr,qliqmsc(1,L),nchp,nchptot,fracland,
13fa5cb63c Jean*1163 & qliqtmp(1,1,L),im,jm)
1662f365b2 Andr*1164 enddo
1165
5317312052 Jean*1166
1167
1662f365b2 Andr*1168 call ctei ( tz,qz,fcctmp,qliqtmp,plz,pkz,pkht,im*jm,nlay )
1169
c909a9788d Andr*1170
5317312052 Jean*1171
1662f365b2 Andr*1172 do l = 1,nlay
5317312052 Jean*1173 do i = 1,ijall
1174 fccavesw(i,1,L) = fccavesw(i,1,L) + fcctmp(i,1,L)
1175 fccavelw(i,1,L) = fccavelw(i,1,L) + fcctmp(i,1,L)
1176 qliqavelw(i,1,L) = qliqavelw(i,1,L) + qliqtmp(i,1,L)
1177 qliqavesw(i,1,L) = qliqavesw(i,1,L) + qliqtmp(i,1,L)
1662f365b2 Andr*1178 enddo
c9f33fa462 Andr*1179 enddo
5317312052 Jean*1180 tmpFac = 1.e6
1181 CALL DIAGNOSTICS_FILL(fcctmp,'TRBFCC ',0,nlay,3,bi,bj,myid)
1182 CALL DIAGNOSTICS_SCALE_FILL( qliqtmp,tmpFac,1,
1183 & 'TRBQLIQ ',0,nlay,3,bi,bj,myid)
1662f365b2 Andr*1184
1185
1186
5317312052 Jean*1187 do i = 1,ijall
1188 tgz(i,1) = 0.
1662f365b2 Andr*1189 enddo
c88fa11306 Andr*1190 call msc2grd(igrd,chfr,tground ,nchp,nchptot,fracland,tgz ,im,jm)
1662f365b2 Andr*1191
5317312052 Jean*1192
1193
1194
1195
1196
13fa5cb63c Jean*1197 call diagnostics_fill(tgz,'TGROUND ',0,1,3,bi,bj,myid)
1198 call diagnostics_fill(tgz,'TCANOPY ',0,1,3,bi,bj,myid)
1662f365b2 Andr*1199
3a3c653dfd Andr*1200 if(diagnostics_is_on('TS ',myid) ) then
5317312052 Jean*1201 do i = 1,ijall
1202 tmpdiag(i,1) = tz(i,1,nlay) * pkht(i,1,nlay)
3a3c653dfd Andr*1203 enddo
1204 call diagnostics_fill(tmpdiag,'TS ',0,1,3,bi,bj,myid)
1662f365b2 Andr*1205 endif
1206
13fa5cb63c Jean*1207 call diagnostics_fill(tgz,'DTG ',0,1,3,bi,bj,myid)
1662f365b2 Andr*1208
5317312052 Jean*1209
1210
1211
1212 tmpFac = atimstp * secday
1213 CALL DIAGNOSTICS_SCALE_FILL(dvturb,tmpFac,1,
1214 & 'TURBV ',0,nlay,3,bi,bj,myid)
1215 CALL DIAGNOSTICS_SCALE_FILL(duturb,tmpFac,1,
1216 & 'TURBU ',0,nlay,3,bi,bj,myid)
1217 tmpFac = atimstp * secday * 1000.
1218 CALL DIAGNOSTICS_SCALE_FILL(dqturb,tmpFac,1,
1219 & 'TURBQ ',0,nlay,3,bi,bj,myid)
1220
13fa5cb63c Jean*1221 if(diagnostics_is_on('TURBT ',myid) ) then
1222 do L = 1,nlay
1223 do i = 1,ijall
5317312052 Jean*1224 tmpdiag(i,1) = dtturb(i,1,l) * pkz(i,1,l)*atimstp*secday
3a3c653dfd Andr*1225 enddo
13fa5cb63c Jean*1226 call diagnostics_fill(tmpdiag,'TURBT ',L,1,3,bi,bj,myid)
1227 enddo
1228 endif
1662f365b2 Andr*1229
5317312052 Jean*1230
1231
1662f365b2 Andr*1232 do i =1,ijall
1233 thtgz(i) = pz(i,1) * atimstp
1234 enddo
1235 do l =1,nlay
1236 do i =1,ijall
1237 duturb(i,1,l) = duturb(i,1,l)*atimstp
1238 dvturb(i,1,l) = dvturb(i,1,l)*atimstp
af35a8ab47 Andr*1239 dtturb(i,1,l) = dtturb(i,1,l)*thtgz(i)
1662f365b2 Andr*1240 enddo
1241 do nt = 1,ntracers
1242 do i =1,ijall
af35a8ab47 Andr*1243 dqturb(i,1,l,nt) = dqturb(i,1,l,nt)*thtgz(i)
1662f365b2 Andr*1244 enddo
1245 enddo
1246 enddo
1247
5317312052 Jean*1248
1249
1250
1662f365b2 Andr*1251
1252 if( time_left.lt.timstp ) then
5317312052 Jean*1253 do i =1,ijall
1254 rainlsp(i,1) = 0.
1255 rainconv(i,1) = 0.
1256 snowfall(i,1) = 0.
1257 enddo
1662f365b2 Andr*1258 endif
1259
1260 return
1261 end
13fa5cb63c Jean*1262
1263
1264
1662f365b2 Andr*1265 SUBROUTINE TRBFLX (NN,TH,THV,SH,U,V,QQ,PL,PLE,PLK,PLKE,DPSTR,
1266 1 IWATER,Z0,tracers,ntrace,ntracedim,DTAU,ITRTRB,KMBG,KHBG,QBEG,
1267 2 TPROF,WU,WV,SRI,ET,EU,SWINDS,sustar,sz0,freqdg,pbldpth,
1268 3 sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
d27a65ec88 Andr*1269 4 irun,nlev,nltop,NYMD,NHMS,grav,cp,rgas,faceps,virtcon,undef,
1662f365b2 Andr*1270 5 dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,
1271 6 hsturb,dhsdqa,dhsdtc,transth,transsh,
1272 7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
5317312052 Jean*1293
1662f365b2 Andr*1294
1295
257b288583 Andr*1296
1297
1662f365b2 Andr*1298
1299
1300
d27a65ec88 Andr*1301
1662f365b2 Andr*1302
1303
1304
1305
1306
1307
1308
175684e43e Andr*1309 implicit none
1662f365b2 Andr*1310
175684e43e Andr*1311
d27a65ec88 Andr*1312 integer nn,irun,nlev,nltop,ntrace,ntracedim,itrtrb,nhms,nymd
a456aa407c Andr*1313 _RL TH(irun,NLEV+1),THV(irun,NLEV+1),SH(irun,NLEV+1)
1314 _RL U(irun,NLEV+1),V(irun,NLEV+1),QQ(irun,NLEV)
1315 _RL PL(irun,NLEV),PLE(irun,NLEV+1),PLK(irun,NLEV)
1316 _RL PLKE(irun,NLEV+1),DPSTR(irun,NLEV)
175684e43e Andr*1317 integer IWATER(irun)
a456aa407c Andr*1318 _RL Z0(irun)
1319 _RL tracers(irun,nlev+1,ntracedim)
1320 _RL dtau,KMBG,KHBG
1662f365b2 Andr*1321 LOGICAL QBEG,TPROF
a456aa407c Andr*1322 _RL SWINDS(irun)
1323 _RL SRI(irun,nlev), ET(irun,nlev)
1324 _RL EU (irun,nlev)
1325 _RL WU(irun,nlev)
1326 _RL WV (irun,nlev), pbldpth(irun)
1327 _RL sustar(irun), sz0(irun)
1328 _RL freqdg(irun,nlev-1)
1329 _RL sct(irun), scu(irun)
1330 _RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
1331 _RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
1332 _RL grav,cp,rgas,faceps,virtcon,undef
1333 _RL eturb(irun),dedqa(irun),dedtc(irun)
1334 _RL hsturb(irun),dhsdqa(irun),dhsdtc(irun)
1335 _RL dshdthg(irun,nlev),dthdthg(irun,nlev)
1336 _RL dshdshg(irun,nlev),dthdshg(irun,nlev)
1337 _RL transth(irun,nlev),transsh(irun,nlev)
1338 _RL ctsave(irun),xxsave(irun),yysave(irun)
1339 _RL zetasave(irun),xlsave(irun,nlev),khsave(irun,nlev)
1340 _RL qliq(irun,nlev),turbfcc(irun,nlev)
1662f365b2 Andr*1341
453e05dab3 Andr*1342
a456aa407c Andr*1343 _RL b1,b3,alpha,halpha,qqmin,qbustr
5317312052 Jean*1344 PARAMETER ( B1 = 16.6 )
1345 PARAMETER ( B3 = 1. / B1 )
175684e43e Andr*1346 PARAMETER ( ALPHA = 0.1 )
1347 PARAMETER ( HALPHA = ALPHA * 0.5 )
1348 PARAMETER ( QQMIN = 0.005 )
1349 PARAMETER ( QBUSTR = 2.550952 )
a456aa407c Andr*1350 _RL argmax, onethrd, z1pem25, b2, two
175684e43e Andr*1351 PARAMETER (ARGMAX = 30.)
1352 PARAMETER (ONETHRD = 1./3. )
1353 PARAMETER (Z1PEM25 = 1.E-25)
1354 PARAMETER ( B2 = 10.1 )
1355 PARAMETER ( two = 2.0 )
1356
a456aa407c Andr*1357 _RL AHS (irun), HS(irun)
1358 _RL XX (irun), YY(irun), CU(irun)
1359 _RL CT(irun), USTAR(irun)
1360 _RL RIB(irun), ZETA(irun), WS(irun)
1361 _RL DTHS(irun), DELTHS(irun)
1362 _RL DTHL(irun), DELTHL(irun)
1363 _RL RIBIN(irun),CUIN(irun)
1364 _RL CTIN(irun),ZETAIN(irun)
1365 _RL USTARIN(irun),RHOSIN(irun),Z0IN(irun)
1366 _RL qqcolmin(irun),qqcolmax(irun),levpbl(irun)
1367
1368 _RL ADZ1(irun,nlev), DZ1TMP(irun,nlev)
1369 _RL DZ3(irun,nlev), TEMP(irun,nlev)
1370 _RL DV(irun,nlev), DTHV(irun,nlev)
1371 _RL DPK(irun,nlev), STRT(irun,nlev)
1372 _RL DW2(irun,nlev), RI(irun,nlev)
1373 _RL RHOZPK(irun,nlev), Q(irun,nlev)
1374 _RL RIINIT(irun,nlev), DU(irun,nlev)
1375 _RL QQINIT(irun,nlev), RHOKDZ(irun,nlev)
1376 _RL RHODZ2(irun,nlev)
1377 _RL KM(irun,nlev), KH(irun,nlev)
1378
1379 _RL DELTH (irun,nlev+1), DELSH (irun,nlev+1)
1380 _RL FLXFAC (irun,nlev+1)
1381 _RL FLXFPK (irun,nlev+1)
1382
1383 _RL ADZ2 (irun,nlev-1), RHODZ1(irun,nlev-1)
1384 _RL VKZE (irun,nlev-1), VKZM (irun,nlev-1)
1385 _RL XL (irun,nlev-1), QXLM (irun,nlev-1)
1386 _RL QQE (irun,nlev-1), QE (irun,nlev-1)
1387 _RL P3 (irun,nlev-1), XQ (irun,nlev-1)
1388 _RL XLDIAG (irun,nlev-1), FLXFCE(irun,nlev-1)
1662f365b2 Andr*1389
1390 LOGICAL FIRST,LAST
453e05dab3 Andr*1391 integer IBITSTB(irun,nlev),INTQ(irun,nlev)
5317312052 Jean*1392
1662f365b2 Andr*1393
1394
a456aa407c Andr*1395 _RL TL(irun,NLEV),DTH(irun,NLEV)
1396 _RL DSH(irun,NLEV)
1397 _RL SHL(irun,NLEV)
1398 _RL AA(irun,NLEV),BB(irun,NLEV),SSDEV(irun,NLEV)
1399 _RL ARG(irun,NLEV),XXZETA(irun),QBYU(irun)
1400 _RL SVAR(irun,NLEV),Q1M(irun,NLEV)
1401 _RL FCC(irun,NLEV)
1402 _RL BETAT(irun,NLEV),BETAW(irun,NLEV)
1403 _RL BETAL(irun,NLEV),BETAT1(irun,NLEV)
1404 _RL BETAW1(irun,NLEV),SBAR(irun,NLEV)
1405 _RL SHSAT(irun,NLEV)
1662f365b2 Andr*1406
1407
1408 logical LWATER
1409 integer IVBITRIB(irun)
a456aa407c Andr*1410 _RL VHZ(irun)
1411 _RL VH0(irun)
1412 _RL VPSIM(irun),VAPSIM(irun)
1413 _RL VPSIG(irun),VPSIHG(irun)
1414 _RL VTEMP(irun),VDZETA(irun)
1415 _RL VDZ0(irun),VDPSIM(irun)
1416 _RL VDPSIH(irun),VZH(irun)
1417 _RL VXX0(irun),VYY0(irun)
1418 _RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
1419 _RL VPSIH(irun),VZETAL(irun)
1420 _RL VZ0L(irun),VPSIH2(irun)
1421 _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
1422 _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
1423 _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
1424 _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
1425 _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
1426 _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
1427 _RL VDPSIMC(irun),VDPSIHC(irun)
1428
1429 _RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
1430 _RL XL0(irun,nlev),Q1(irun,nlev-1)
1431 _RL WRKIT1(irun,nlev-1)
1432 _RL WRKIT2(irun,nlev-1)
1433 _RL WRKIT3(irun,nlev-1)
1434 _RL WRKIT4(irun,nlev-1)
1662f365b2 Andr*1435 INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
5317312052 Jean*1436
a456aa407c Andr*1437 _RL vrt1con,pi,rsq2pi,p5sr,clh,vk,rvk,aitr,gbycp,fac1,fac2
1438 _RL getcon,dum,errf
453e05dab3 Andr*1439 integer istnlv,nlevm1,nlevm2,nlevml,nlevp1,istnm1,istnm2,istnp1
1440 integer istnml,istnmq,istlmq,nlevmq
45ad04df92 Jean*1441 integer i,iter,init,n,LL,L,Lp,Lp1,lmin,lminq,lminq1,ibit
1662f365b2 Andr*1442
1443 vk = getcon('VON KARMAN')
1444 rvk = 1./vk
1445 AITR = 1. / FLOAT(ITRTRB)
1446 ISTNLV = irun * NLEV
1447 NLEVM1 = NLEV - 1
1448 NLEVM2 = NLEV - 2
1449 NLEVP1 = NLEV + 1
1450 ISTNM1 = irun * NLEVM1
1451 ISTNM2 = irun * NLEVM2
1452 ISTNP1 = irun * NLEVP1
1453 GBYCP = GRAV / CP
5317312052 Jean*1454
1662f365b2 Andr*1455 VRT1CON = 1. + VIRTCON
1456 PI = 4. * ATAN(1.)
1457 RSQ2PI = 1./ ((2.*PI)**0.5)
1458 P5SR = 0.5**0.5
1459 CLH = GETCON('LATENT HEAT COND') / CP
1460
1461
1462
1463 N = 6
1464
1465
1466 INIT = 0
1467 IF(QBEG) INIT = 1
1468
1469
5317312052 Jean*1470 DO L =1,NLEV
1471 DO I =1,irun
1472 wu(I,L) = 0.
1473 wv(I,L) = 0.
1474 eu(I,L) = 0.
1475 et(I,L) = 0.
1476 ENDDO
1477 ENDDO
1662f365b2 Andr*1478 if (tprof) then
5317312052 Jean*1479 DO L =1,NLEVM1
1480 DO I =1,irun
1481 XLDIAG(I,L) = 0.
1482 ENDDO
1483 ENDDO
1662f365b2 Andr*1484 endif
5317312052 Jean*1485 DO L =1,NLEVM1
1486 DO I =1,irun
1487 FREQDG(I,L) = 0.
1488 ENDDO
1489 ENDDO
1662f365b2 Andr*1490 do I =1,irun
1491 scu(i) = 0.
1492 enddo
1493 do I =1,irun
1494 sct(i) = 0.
1495 enddo
1496 do I =1,irun
1497 pbldpth(i) = 0.
1498 enddo
1499 do I =1,irun
1500 sustar(i) = 0.
1501 enddo
1502 do I =1,irun
1503 sz0(i) = 0.
1504 enddo
5317312052 Jean*1505
1506
1507
1662f365b2 Andr*1508 do I =1,irun
1509 stu2m(i) = 0.
1510 enddo
1511 do I =1,irun
1512 stv2m(i) = 0.
1513 enddo
1514 do I =1,irun
1515 stt2m(i) = 0.
1516 enddo
1517 do I =1,irun
1518 stq2m(i) = 0.
1519 enddo
1520 do I =1,irun
1521 stu10m(i) = 0.
1522 enddo
1523 do I =1,irun
1524 stv10m(i) = 0.
1525 enddo
1526 do I =1,irun
1527 stt10m(i) = 0.
1528 enddo
1529 do I =1,irun
1530 stq10m(i) = 0.
1531 enddo
1532
1533 IF (INIT.EQ.1) THEN
5317312052 Jean*1534 DO L =1,NLEVM1
1535 DO I =1,irun
1536 XLSAVE(I,L) = 0.
1537 KHSAVE(I,L) = 0.
1538 ENDDO
1662f365b2 Andr*1539 ENDDO
1540 DO I = 1,irun
1541 CTSAVE(I) = 0.
1542 XXSAVE(I) = 0.
1543 YYSAVE(I) = 0.
1544 ZETASAVE(I) = 0.
1545 ENDDO
1546 ENDIF
1547
1548
1549
5317312052 Jean*1550 DO L =1,NLEV
1551 DO I =1,irun
1552 ADZ1(I,L) = (CP/GRAV)*(PLKE(I,L+1)-PLKE(I,L))
1553 ADZ1(I,L) = THV(I,L) * ADZ1(I,L)
1554 DZ1TMP(I,L) = ADZ1(I,L)
1555 ENDDO
1556 ENDDO
1557 DO L =1,NLEVM1
1558 DO I =1,irun
1559 ADZ2(I,L) = 0.5 * (ADZ1(I,L)+ADZ1(I,L+1))
1560 ENDDO
1561 ENDDO
1662f365b2 Andr*1562
1563
1564 DO 9042 I =1,irun
1565 HS(I) = 0.5 * ADZ1(I,NLEV)
1566 9042 CONTINUE
1567
1568
1569 DO 9044 I =1,irun
1570 DZ3(I,1) = HALPHA * ADZ1(I,1)
1571 9044 CONTINUE
5317312052 Jean*1572 DO L =2,NLEVM1
1573 DO I =1,irun
1574 DZ3(I,L) = ALPHA * ADZ1(I,L)
1575 ENDDO
1576 ENDDO
1662f365b2 Andr*1577 DO 9048 I =1,irun
1578 DZ3(I,NLEV) = ALPHA * HS(I)
1579 9048 CONTINUE
1580
1581
1582
5317312052 Jean*1583 DO L =2,NLEV
1584 DO I =1,irun
1585 TEMP(I,L) = VK * ADZ1(I,L)
1586 ENDDO
1587 ENDDO
1662f365b2 Andr*1588 DO 9052 I =1,irun
1589 VKZE(I,NLEVM1) = TEMP(I,NLEV)
1590 9052 CONTINUE
1591 DO 100 LL = 2,NLEVM1
1592 L = NLEV - LL
1593 LP1 = L + 1
1594 DO 9054 I =1,irun
1595 VKZE(I,L) = VKZE(I,LP1) + TEMP(I,LP1)
1596 9054 CONTINUE
1597 100 CONTINUE
5317312052 Jean*1598 DO L =1,NLEVM1
1599 DO I =1,irun
1600 VKZM(I,L) = VKZE(I,L) - 0.5 * TEMP(I,L+1)
1601 ENDDO
1602 ENDDO
1662f365b2 Andr*1603
1604
1605 DO 200 L = 1,NLEVM1
1606 LP1 = L + 1
1607 DO 9058 I =1,irun
1608 FAC1 = DPSTR(I,L) / ( DPSTR(I,L) + DPSTR(I,LP1) )
1609 FAC2 = 1. - FAC1
1610 RHODZ2(I,L) = FAC1 * THV(I,LP1)
1611 RHODZ2(I,L) = RHODZ2(I,L) + FAC2 * THV(I,L)
1612 9058 CONTINUE
1613 200 CONTINUE
5317312052 Jean*1614 DO L =1,NLEVM1
1615 DO I =1,irun
1616 RHODZ2(I,L) = (RGAS*0.01) * RHODZ2(I,L)
1617 TEMP(I,L) = PLKE(I,L+1) * ADZ2(I,L)
1618 RHODZ2(I,L) = TEMP(I,L) * RHODZ2(I,L)
1619 RHODZ2(I,L) = PLE(I,L+1) / RHODZ2(I,L)
1620 RHOZPK(I,L) = RHODZ2(I,L) * PLKE(I,L+1)
1621 RHODZ1(I,L) = (RGAS*0.01) * THV(I,L+1)
1622 TEMP(I,L) = PLK(I,L+1) * ADZ1(I,L+1)
1623 RHODZ1(I,L) = TEMP(I,L) * RHODZ1(I,L)
1624 RHODZ1(I,L) = PL(I,L+1) / RHODZ1(I,L)
1625 ENDDO
1626 ENDDO
1662f365b2 Andr*1627
1628
1629
5317312052 Jean*1630 DO L =1,NLEV
1631 DO I =1,irun
1632 FLXFPK(I,L) = PLE(I,L+1) - PLE(I,L)
1633 FLXFPK(I,L) = FLXFPK(I,L) * PLK(I,L)
1634 FLXFPK(I,L) = (GRAV*DTAU*0.01) / FLXFPK(I,L)
1635 ENDDO
1636 ENDDO
1662f365b2 Andr*1637 DO 9064 I =1,irun
1638 FLXFPK(I,NLEVP1) = 0.
1639 9064 CONTINUE
1640 DO 9066 I =1,irun
1641 IF (IWATER(I).EQ.0 ) FLXFPK(I,NLEVP1) = 1. / PLKE(I,NLEVP1)
1642 9066 CONTINUE
5317312052 Jean*1643 DO L =1,NLEV
1644 DO I =1,irun
1645 FLXFAC(I,L) = FLXFPK(I,L) * PLK(I,L)
1646 ENDDO
1647 ENDDO
1662f365b2 Andr*1648 DO 9070 I =1,irun
1649 FLXFAC(I,NLEVP1) = FLXFPK(I,NLEVP1)
1650 9070 CONTINUE
1651 DO 9074 I =1,irun
1652 FLXFPK(I,NLEVP1) = CP * FLXFPK(I,NLEVP1)
1653 9074 CONTINUE
5317312052 Jean*1654 DO L =1,NLEVM1
1655 DO I =1,irun
1656 FLXFCE(I,L) = PL(I,L+1) - PL(I,L)
1657 FLXFCE(I,L) = (GRAV*DTAU*0.01) / FLXFCE(I,L)
1658 ENDDO
1659 ENDDO
1662f365b2 Andr*1660
1661
5317312052 Jean*1662 DO L =1,NLEV
1663 DO I =1,irun
1664 ADZ1(I,L) = 1. / ADZ1(I,L)
1665 ENDDO
1666 ENDDO
1667 DO L =1,NLEVM1
1668 DO I =1,irun
1669 ADZ2(I,L) = 1. / ADZ2(I,L)
1670 ENDDO
1671 ENDDO
1662f365b2 Andr*1672 DO 9088 I =1,irun
1673 AHS(I) = 1. / HS(I)
1674 9088 CONTINUE
1675
1676
5317312052 Jean*1677 DO L =1,NLEVM1
1678 DO I =1,irun
1679 DPK(I,L) = ( PLK(I,L+1)-PLK(I,L) ) * ADZ2(I,L)
1680 ENDDO
1681 ENDDO
1662f365b2 Andr*1682 DO 9092 I =1,irun
1683 DPK(I,NLEV) = GBYCP / THV(I,NLEV)
1684 9092 CONTINUE
1685
1686
5317312052 Jean*1687 DO L =1,NLEVM1
1688 DO I =1,irun
1689 Q(I,L) = 2. * QQ(I,L)
1690 Q(I,L) = SQRT( Q(I,L) )
1691 ENDDO
1692 ENDDO
1662f365b2 Andr*1693 FIRST = .TRUE.
1694 LAST = .FALSE.
1695
1696
1697
13fa5cb63c Jean*1698
1662f365b2 Andr*1699 DO 2000 ITER = 1, ITRTRB
13fa5cb63c Jean*1700
1662f365b2 Andr*1701 IF ( ITER .GE. ITRTRB ) LAST = .TRUE.
13fa5cb63c Jean*1702
1662f365b2 Andr*1703
13fa5cb63c Jean*1704
1662f365b2 Andr*1705 IF(ITER.EQ.1) THEN
1706 DO I = 1,irun
1707 CT(I) = CTSAVE(I)
1708 XX(I) = XXSAVE(I)
1709 YY(I) = YYSAVE(I)
1710 ZETA(I) = ZETASAVE(I)
1711 ENDDO
1712 ENDIF
13fa5cb63c Jean*1713
1662f365b2 Andr*1714 DO I = 1,irun
1715 TL(I,NLEV) = TH(I,NLEV)*PLK(I,NLEV)
1716 call qsat ( tl(i,nlev),pl(i,nlev),shsat(i,nlev),dum,.false. )
1717 ENDDO
1718
1719 DO I = 1,irun
1720 BB(I,NLEV) = FACEPS*SHSAT(I,NLEV)/(TL(I,NLEV)*TL(I,NLEV))
1721 AA(I,NLEV) = 1. / (1. + CLH * BB(I,NLEV) )
1722 BB(I,NLEV) = BB(I,NLEV) * AA(I,NLEV) * plk(I,nlev)
1723 DTH(I,NLEV) = TH(I,NLEV)-TH(I,NLEVP1)
1724 DSH(I,NLEV) = SH(I,NLEV)-SH(I,NLEVP1)
1725 SBAR(I,NLEV) = AA(I,NLEV) * (SH(I,NLEV) - SHSAT(I,NLEV))
1726 SSDEV(I,NLEV)=CT(I)*(AA(I,NLEV)*DSH(I,NLEV)
1727 1 -BB(I,NLEV)*DTH(I,NLEV))
1728 XXZETA(I) = XX(I)-ZETA(I)
1729 IF(XXZETA(I).LT.0.1*XX(I)) XXZETA(I)=0.1*XX(I)
1730 IF(XXZETA(I).LE.0.) XXZETA(I)=0.1
1731 QBYU(I) =QBUSTR * XXZETA(I) ** ONETHRD
1732 SSDEV(I,NLEV) = B2*YY(I)*SSDEV(I,NLEV)*SSDEV(I,NLEV)/QBYU(I)
1733 SVAR(I,NLEV) = SQRT(SSDEV(I,NLEV))
1734 IF ( SVAR(I,NLEV).LT.Z1PEM25) SVAR(I,NLEV) = Z1PEM25
1735 Q1M(I,NLEV) = SBAR(I,NLEV) / SVAR(I,NLEV)
1736 FCC(I,NLEV) = (1./2.) * ( 1. + ERRF( P5SR*Q1M(I,NLEV) ) )
1737 SHL(I,NLEV) = FCC(I,NLEV) * SBAR(I,NLEV)
1738 ARG(I,NLEV) = (1./2.)*Q1M(I,NLEV)*Q1M(I,NLEV)
1739 IF(ARG(I,NLEV).LE.ARGMAX)
1740 1 SHL(I,NLEV) = SHL(I,NLEV)+RSQ2PI*SVAR(I,NLEV)*EXP(-ARG(I,NLEV))
1741 BETAT(I,NLEV) = 1. + VIRTCON*SH(I,NLEV) - VRT1CON*SHL(I,NLEV)
1742 BETAW(I,NLEV) = VIRTCON *
1743 1 ( TH(I,NLEV) + CLH * SHL(I,NLEV) * (1./plk(i,nlev)) )
1744 BETAL(I,NLEV) = (1.+VIRTCON*SH(I,NLEV)-TWO*VRT1CON*SHL(I,NLEV))
1745 1 * (1./plk(i,nlev)) * CLH - VRT1CON * TH(I,NLEV)
1746 BETAT1(I,NLEV) = BETAT(I,NLEV) - BB(I,NLEV)*FCC(I,NLEV)
1747 1 * BETAL(I,NLEV)
1748 BETAW1(I,NLEV) = BETAW(I,NLEV) + AA(I,NLEV) * FCC(I,NLEV)
1749 1 * BETAL(I,NLEV)
1750 DTHV(I,NLEV) = BETAT1(I,NLEV)*DTH(I,NLEV) +
1751 1 BETAW1(I,NLEV)*DSH(I,NLEV)
1752 THV(I,NLEVP1) = THV(I,NLEV) - DTHV(I,NLEV)
1753 ENDDO
1754
1755
13fa5cb63c Jean*1756
1662f365b2 Andr*1757 CALL SFCFLX(NN,U(1,NLEV),V(1,NLEV),
1758 1 THV(1,NLEV),
1759 2 THV(1,NLEVP1),TH(1,NLEV),TH(1,NLEVP1),
1760 3 SH(1,NLEV),SH(1,NLEVP1),PLK(1,NLEV),
1761 4 PLKE(1,NLEVP1),PLE(1,NLEVP1),Z0,
1762 5 IWATER,HS,AHS,
1763 6 FIRST,LAST,N,irun,aitr,RHODZ2(1,NLEV),RHOZPK(1,NLEV),
1764 7 KH(1,NLEV),KM(1,NLEV),USTAR,
1765 8 XX,YY,CU,
1766 9 CT,RIB,ZETA,WS,
1767 1 stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
1768 2 cp,rgas,undef,
1769 3 lwater, ivbitrib,
1770 4 VHZ,VPSIM,VAPSIM,VPSIG,VPSIHG,VTEMP,VDZETA,VDZ0,VDPSIM,
1771 5 VDPSIH,VZH,VXX0,VYY0,VAPSIHG,VRIB1,VWS1,VPSIH,
1772 9 VZETAL,VZ0L,VPSIH2,VH0,
1773 1 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
1774 2 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
1775 3 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
1776
13fa5cb63c Jean*1777
1662f365b2 Andr*1778 N = 1
13fa5cb63c Jean*1779
1662f365b2 Andr*1780
13fa5cb63c Jean*1781
1662f365b2 Andr*1782
1783 DO 9098 I =1,irun
1784 Q(I,NLEV) = QBUSTR * USTAR(I)
1785 QQ(I,NLEV) = 0.5 * Q(I,NLEV) * Q(I,NLEV)
1786 9098 CONTINUE
1787
13fa5cb63c Jean*1788
1662f365b2 Andr*1789
1790
5317312052 Jean*1791 DO L =1,NLEVM1
1792 DO I =1,irun
1793 DU(I,L) = ( U(I,L)- U(I,L+1) ) * ADZ2(I,L)
1794 DV(I,L) = ( V(I,L)- V(I,L+1) ) * ADZ2(I,L)
1795 ENDDO
1796 ENDDO
1662f365b2 Andr*1797
1798
13fa5cb63c Jean*1799
5317312052 Jean*1800 IF(ITER.EQ.1) THEN
1801 DO L =1,NLEVM1
1802 DO I =1,irun
1803 XL(I,L) = XLSAVE(I,L)
1804 ENDDO
1662f365b2 Andr*1805 ENDDO
1806 ENDIF
13fa5cb63c Jean*1807
5317312052 Jean*1808 DO L =1,NLEVM1
1809 DO I =1,irun
1810 DTH(I,L) = ( TH(I,L)-TH(I,L+1) ) * ADZ2(I,L)
1811 DSH(I,L) = ( SH(I,L)-SH(I,L+1) ) * ADZ2(I,L)
1812 TL(I,L) = TH(I,L)*PLK(I,L)
1813 ENDDO
1662f365b2 Andr*1814 ENDDO
5317312052 Jean*1815 DO LL = 1,NLEVM1
1816 DO I = 1,irun
1662f365b2 Andr*1817 call qsat ( tl(i,LL),pl(i,LL),shsat(i,LL),dum,.false. )
1818 ENDDO
1819 ENDDO
5317312052 Jean*1820 DO L =1,NLEVM1
1821 DO I =1,irun
1822 BB(I,L) = FACEPS*SHSAT(I,L)/(TL(I,L)*TL(I,L))
1823 AA(I,L) = 1. / (1. + CLH * BB(I,L) )
1824
1825 BB(I,L) = BB(I,L) * AA(I,L)
1826 SBAR(I,L) = AA(I,L) * (SH(I,L) - SHSAT(I,L))
1827 ENDDO
1662f365b2 Andr*1828 ENDDO
1829 DO I = 1,irun
1830
1831 SSDEV(I,1) = XL(I,1)*(AA(I,1)*DSH(I,1)-
1832 1 BB(I,1)*plke(I,2)*DTH(I,1))
1833 SSDEV(I,1) = B2 * KHSAVE(I,1) * SSDEV(I,1) * SSDEV(I,1)
1834 SVAR(I,1) = SQRT(SSDEV(I,1))
1835 IF ( SVAR(I,1).LT.Z1PEM25) SVAR(I,1) = Z1PEM25
1836 ENDDO
5317312052 Jean*1837 DO L =2,NLEVM1
1838 DO I =1,irun
1839
1840 SSDEV(I,L) = XL(I,L-1)*(AA(I,L)*DSH(I,L-1)-
1841 1 BB(I,L)*plke(I,L)*DTH(I,L-1))
1842 SSDEV(I,L) = B2 * KHSAVE(I,L-1) * SSDEV(I,L-1) * SSDEV(I,L-1)
1843 SVAR(I,L) = SQRT(SSDEV(I,L))
1844
1845 SSDEV(I,L) = XL(I,L)*(AA(I,L)*DSH(I,L)-
1846 1 BB(I,L)*plke(I,L+1)*DTH(I,L))
1847 SSDEV(I,L) = B2 * KHSAVE(I,L) * SSDEV(I,L) * SSDEV(I,L)
1848 TEMP(I,L) = SQRT(SSDEV(I,L))
1849 SVAR(I,L) = (1./2.) * (SVAR(I,L) + TEMP(I,L))
1850 IF ( SVAR(I,L).LT.Z1PEM25) SVAR(I,L) = Z1PEM25
1851 ENDDO
1662f365b2 Andr*1852 ENDDO
5317312052 Jean*1853 DO L =1,NLEVM1
1854 DO I =1,irun
1855 Q1M(I,L) = SBAR(I,L) / SVAR(I,L)
1856 FCC(I,L) = (1./2.) * ( 1. + ERRF( P5SR*Q1M(I,L) ) )
1857 SHL(I,L) = FCC(I,L) * SBAR(I,L)
1858 ARG(I,L) = (1./2.)*Q1M(I,L)*Q1M(I,L)
1859 IF(ARG(I,L).LE.ARGMAX)
1860 1 SHL(I,L) = SHL(I,L)+RSQ2PI*SVAR(I,L)*EXP(-ARG(I,L))
1861 BETAT(I,L) = 1. + VIRTCON * SH(I,L) - VRT1CON * SHL(I,L)
1862 BETAW(I,L) = VIRTCON *
1863 1 ( TH(I,L) + (CLH/plk(I,L)) * SHL(I,L) )
1864 BETAL(I,L) = ( 1. + VIRTCON*SH(I,L) - TWO*VRT1CON*SHL(I,L) )
1865 1 * (CLH/plke(I,L+1)) - VRT1CON * TH(I,L)
1866
1867 BETAT1(I,L) = BETAT(I,L) -
1868 1 BB(I,L)*plk(i,L) * FCC(I,L) * BETAL(I,L)
1869 BETAW1(I,L) = BETAW(I,L) + AA(I,L) * FCC(I,L) * BETAL(I,L)
1870 ENDDO
1662f365b2 Andr*1871 ENDDO
5317312052 Jean*1872 DO L =1,NLEVM1
1873 DO I =1,irun
1874 DTHV(I,L) = (1./2.)*((BETAT1(I,L)+BETAT1(I,L+1))*DTH(I,L)
1875 1 + (BETAW1(I,L)+BETAW1(I,L+1))*DSH(I,L))
1876 ENDDO
1662f365b2 Andr*1877 ENDDO
5317312052 Jean*1878
1662f365b2 Andr*1879
1880
1881 DO 9102 I =1,irun
1882 DU(I,NLEV) = CU(I)*XX(I)*AHS(I)*RVK
1883 DV(I,NLEV) = V(I,NLEV) * DU(I,NLEV)
1884 DU(I,NLEV) = U(I,NLEV) * DU(I,NLEV)
1885 DTHV(I,NLEV) = CT(I) * YY(I) *
1886 1 ((THV(I,NLEV)-THV(I,NLEVP1)) * RVK)* AHS(I)
1887 9102 CONTINUE
1888
1889
1890
5317312052 Jean*1891 DO L =1,NLEV
1892 DO I =1,irun
1893 STRT(I,L) = CP * DTHV(I,L) * DPK(I,L)
1894 DW2(I,L) = DU(I,L) * DU(I,L) + DV(I,L) * DV(I,L)
1895 IF ( DW2(I,L) .LE. 1.e-4 ) DW2(I,L) = 1.e-4
1896 RI(I,L) = STRT(I,L) / DW2(I,L)
1897 ENDDO
1898 ENDDO
1662f365b2 Andr*1899
1900
1901
5317312052 Jean*1902 DO L =1,NLEVM1
1903 DO I =1,irun
1904 SRI(I,L) = RI(I,L)
1905 ENDDO
1906 ENDDO
1662f365b2 Andr*1907 DO 9108 I =1,irun
1908 SRI(I,NLEV) = RIB(I)
1909 9108 CONTINUE
1910 DO 9110 I =1,irun
1911 SWINDS(I) = WS(I)
1912 9110 CONTINUE
1913
1914
5317312052 Jean*1915 DO L =1,NLEVM1
1916 DO I =1,irun
1917 KH(I,L) = 0.
1918 KM(I,L) = 0.
1919 QQE(I,L) = 0.
1920 QE(I,L) = 0.
1921 P3(I,L) = 0.
1922 IBITSTB(I,L) = 0
1923 IF ( QQ(I,L) .GT. 1.e-8 ) THEN
1924 INTQ(I,L) = 1
1925 ELSE
1926 INTQ(I,L) = 0
1927 ENDIF
1928 IF ( QQ(I,L).LE.1.e-8 ) THEN
1929 QQ(I,L) = 0.
1930 Q(I,L) = 0.
1931 ENDIF
1932 ENDDO
1933 ENDDO
13fa5cb63c Jean*1934
1662f365b2 Andr*1935 DO 300 LMINQ = 1,NLEVM1
1936 IBIT = 0
1937 DO 9116 I = 1,irun
1938 IF ( QQ(I,LMINQ).GT.1.e-8 ) IBIT = IBIT + 1
1939 9116 CONTINUE
1940 IF(IBIT.GE.1)GO TO 310
1941 300 CONTINUE
1942 LMINQ = NLEV-1
1943 310 CONTINUE
1944 LMINQ = 1
1945 LMINQ1 = 1
1946 IF(LMINQ.GT.1)LMINQ1 = LMINQ - 1
1947
1948
1949 CALL TRBLEN(STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,QXLM,
1950 1 NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
1951 2 DZITRP,STBFCN,XL0,Q1,WRKIT1,WRKIT2,WRKIT3,WRKIT4,irun)
1952
1953
1954 IF( LMIN .LT. NLEV ) THEN
1955 NLEVML = NLEV - LMIN
1956 CALL TRBL20(RI(1,LMIN),STRT(1,LMIN),DW2(1,LMIN),XL(1,LMIN),
1957 1 KM(1,LMIN),KH(1,LMIN),QE(1,LMIN),QQE(1,LMIN),IBITSTB(1,LMIN),
1958 2 NLEVML,nlev,irun)
1959 ENDIF
1960
1961
1962 IF ( INIT .EQ. 1 ) THEN
5317312052 Jean*1963 DO L =1,NLEVM1
1964 DO I =1,irun
1965 QQ(I,L) = QQE(I,L)
1966 Q(I,L) = QE(I,L)
1967 ENDDO
1968 ENDDO
1662f365b2 Andr*1969 INIT = 2
1970 CALL TRBLEN(STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,QXLM,
1971 1 NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
1972 2 DZITRP,STBFCN,XL0,Q1,WRKIT1,WRKIT2,WRKIT3,WRKIT4,irun)
1973 INIT = 0
1974 GO TO 550
1975 ENDIF
1976
1977
1978 IF( LMIN .LT. NLEV ) THEN
1979 ISTNML = irun * NLEVML
5317312052 Jean*1980 DO L =LMIN,NLEVM1
1981 DO I =1,irun
1982 IF ( (IBITSTB(I,L).EQ.1) .AND.
1983 1 ( Q(I,L) .LE. QE(I,L) ) ) THEN
1984 IBITSTB(I,L) = 1
1985 ELSE
1986 IBITSTB(I,L) = 0
1987 ENDIF
1988 ENDDO
1989 ENDDO
1990 DO L =LMIN,NLEVM1
1991 DO I =1,irun
1992 IF(IBITSTB(I,L).EQ.1 ) THEN
1993 TEMP(I,L) = Q(I,L) / QE(I,L)
1994 KH(I,L) = TEMP(I,L) * KH(I,L)
1995 KM(I,L) = TEMP(I,L) * KM(I,L)
1996 ENDIF
1997 TEMP(I,L) = 0.01 * QQE(I,L)
1998 IF((IBITSTB(I,L).EQ.1) .AND.
1999 1 ( QQ(I,L) .LE. TEMP(I,L) )) THEN
2000 QQ(I,L) = TEMP(I,L)
2001 Q(I,L) = 0.1 * QE(I,L)
2002 ENDIF
2003 IF(IBITSTB(I,L).EQ.1 ) P3(I,L) = (2.*B3) *
2004 1 ( QE(I,L) - Q(I,L) )
2005 ENDDO
2006 ENDDO
1662f365b2 Andr*2007 ENDIF
2008
2009
2010 NLEVML = NLEV - LMINQ
2011 CALL TRBL25(Q(1,LMINQ),XL(1,LMINQ),STRT(1,LMINQ),DW2(1,LMINQ),
2012 1 IBITSTB(1,LMINQ),INTQ(1,LMINQ),KM(1,LMINQ),KH(1,LMINQ),
2013 2 P3(1,LMINQ),NLEVML,nlev,irun)
2014
2015
2016 IF ( LMINQ .LT. LMIN ) THEN
2017 LMIN = LMINQ
2018 ISTNML = irun * ( NLEV - LMIN )
2019 ENDIF
2020 IF( LMIN .LT. NLEV ) THEN
5317312052 Jean*2021 DO L =LMIN,NLEVM1
2022 DO I =1,irun
2023 P3(I,L) = P3(I,L) * DTAU / XL(I,L)
2024 TEMP(I,L) = QQE(I,L) * P3(I,L)
2025 XQ(I,L) = QQE(I,L) - QQ(I,L)
2026 ENDDO
2027 ENDDO
2028 DO L =LMIN,NLEVM1
2029 DO I =1,irun
2030 IF( ( (IBITSTB(I,L).EQ.1) .AND.
2031 1 ( XQ(I,L) .LT. TEMP(I,L) ) )
1662f365b2 Andr*2032 2 .OR.
5317312052 Jean*2033 3 ( (IBITSTB(I,L).EQ.0) .AND.
2034 4 ( XQ(I,L) .GT. TEMP(I,L) ) ) )
2035 5 P3(I,L) = XQ(I,L) / QQE(I,L)
2036 ENDDO
2037 ENDDO
1662f365b2 Andr*2038 ENDIF
2039 550 CONTINUE
2040
2041
2042 IF ( TPROF .AND. FIRST ) THEN
2043 DO 9118 I =1,irun
2044 RIBIN(I) = RIB(I)
2045 CUIN(I) = CU(I)
2046 CTIN(I) = CT(I)
2047 USTARIN(I) = USTAR(I)
2048 RHOSIN(I) = RHODZ2(I,NLEV)
2049 Z0IN(I) = Z0(I)
2050 ZETAIN(I) = ZETA(I)
2051 9118 CONTINUE
5317312052 Jean*2052 DO L =1,NLEV
2053 DO I =1,irun
2054 RIINIT(I,L) = RI(I,L)
2055 QQINIT(I,L) = QQ(I,L)
2056 ENDDO
2057 ENDDO
1662f365b2 Andr*2058 ENDIF
5f99f321da Andr*2059
2060
2061 do L = 1,nlev-1
2062 do i = 1,irun
13fa5cb63c Jean*2063
2064
2065
5f99f321da Andr*2066 if(pl(i,L).le.150.) then
010bc8d74c Mart*2067 qxlm(i,L) = min(qxlm(i,L),5. _d 0)
5f99f321da Andr*2068 endif
2069 enddo
2070 enddo
13fa5cb63c Jean*2071
1662f365b2 Andr*2072
2073
2074 NLEVMQ = NLEV - LMINQ1
2075 ISTNMQ = irun * NLEVMQ
5317312052 Jean*2076 DO L =LMINQ1,NLEVM1
2077 DO I =1,irun
2078 RHOKDZ(I,L) = RHODZ1(I,L) * QXLM(I,L)
2079 ENDDO
2080 ENDDO
1662f365b2 Andr*2081 CALL TRBDIF(QQ(1,LMINQ1),P3(1,LMINQ1),RHOKDZ(1,LMINQ1),
daaf4a8d7f Andr*2082 1 FLXFCE(1,LMINQ1),DTHS,DELTHS,NLEVMQ,1,1.0 _d -20,irun)
13fa5cb63c Jean*2083
1662f365b2 Andr*2084
13fa5cb63c Jean*2085
5317312052 Jean*2086 DO L =1,NLEVM1
2087 DO I =1,irun
2088 KHSAVE(I,L) = KH(I,L)
2089 ENDDO
1662f365b2 Andr*2090 ENDDO
13fa5cb63c Jean*2091
1662f365b2 Andr*2092
13fa5cb63c Jean*2093
1662f365b2 Andr*2094 IF(LMINQ1.GT.1)THEN
2095 ISTLMQ = irun * (LMINQ1-1)
2096
5317312052 Jean*2097 DO L =1,LMINQ1-1
2098 DO I =1,irun
2099 KM(I,L) = KMBG
2100 KH(I,L) = KHBG
2101 ENDDO
2102 ENDDO
1662f365b2 Andr*2103
2104 ENDIF
13fa5cb63c Jean*2105
1662f365b2 Andr*2106
5317312052 Jean*2107 DO L =LMINQ1,NLEVM1
2108 DO I =1,irun
2109 Q(I,L) = 2. * QQ(I,L)
2110 Q(I,L) = SQRT(Q(I,L))
2111 XQ(I,L)= XL(I,L) * Q(I,L)
2112 KM(I,L)=XQ(I,L)*KM(I,L)+KMBG
2113 KH(I,L)=XQ(I,L)*KH(I,L)+KHBG
2114 ENDDO
2115 ENDDO
1662f365b2 Andr*2116
5f99f321da Andr*2117
67539bd6f4 Andr*2118 do L = 1,nlev-1
2119 do i = 1,irun
13fa5cb63c Jean*2120
2121
2122
2123
67539bd6f4 Andr*2124 if(pl(i,L).le.150.) then
010bc8d74c Mart*2125 kh(i,L) = min(kh(i,L),5. _d 0)
2126 km(i,L) = min(km(i,L),5. _d 0)
67539bd6f4 Andr*2127 endif
2128 enddo
2129 enddo
13fa5cb63c Jean*2130
1662f365b2 Andr*2131
13fa5cb63c Jean*2132
5317312052 Jean*2133 DO L =1,NLEV
2134 DO I =1,irun
2135 TEMP(I,L) = RHOZPK(I,L) * KH(I,L)
2136 DELTH(I,L) = 0.
2137 ENDDO
2138 ENDDO
1662f365b2 Andr*2139 DO 9132 I =1,irun
2140 DELTH(I,NLEVP1) = 1.
2141 9132 CONTINUE
daaf4a8d7f Andr*2142 CALL TRBDIF(TH,DELTH,TEMP,FLXFPK,DTHS,DELTHS,NLEV,2,0. _d 0,irun)
1662f365b2 Andr*2143 do i = 1,irun
2144 hsturb(i) = -1.* dths(i)
2145 dhsdtc(i) = -1.* delths(i)
2146 enddo
2147 do L = 1,nlev
2148 do i = 1,irun
2149 dthdthg(i,L) = delth(i,L)
2150 enddo
2151 enddo
2152 do L = 1,nlev
2153 do i = 1,irun
2154 transth(i,L) = temp(i,L)
2155 enddo
2156 enddo
2157
5317312052 Jean*2158 DO L =1,NLEV
2159 DO I =1,irun
2160 RHOKDZ(I,L) = RHODZ2(I,L) * KH(I,L)
2161 DELSH(I,L) = 0.
2162 ENDDO
2163 ENDDO
1662f365b2 Andr*2164 DO 9140 I =1,irun
2165 DELSH(I,NLEVP1) = 1.
2166 9140 CONTINUE
2167
daaf4a8d7f Andr*2168 CALL TRBDIF(SH,DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,NLEV,
13fa5cb63c Jean*2169 & 2,0. _d 0,irun)
1662f365b2 Andr*2170 do i = 1,irun
2171 eturb(i) = -1.* dthl(i)
2172 dedqa(i) = -1.* delthl(i)
2173 enddo
2174 do L = 1,nlev
2175 do i = 1,irun
2176 dshdshg(i,L) = delsh(i,L)
2177 enddo
2178 enddo
2179 do L = 1,nlev
2180 do i = 1,irun
2181 transsh(i,L) = rhokdz(i,L)
2182 enddo
2183 enddo
2184
2185
13fa5cb63c Jean*2186
1662f365b2 Andr*2187 do i = 1,irun
5317312052 Jean*2188 rhokdz(i,nlev) = 0.0
1662f365b2 Andr*2189 enddo
2190
257b288583 Andr*2191
2192
2193
2194
2195
2196
2197
13fa5cb63c Jean*2198
1662f365b2 Andr*2199
13fa5cb63c Jean*2200
5317312052 Jean*2201 DO L =1,NLEV
2202 DO I =1,irun
2203 RHOKDZ(I,L) = RHODZ2(I,L) * KM(I,L)
2204 ENDDO
2205 ENDDO
daaf4a8d7f Andr*2206 CALL TRBDIF(U,V,RHOKDZ,FLXFAC,DTHS,DELTHS,NLEV,3,0. _d 0,irun)
1662f365b2 Andr*2207
5317312052 Jean*2208 DO L =1,NLEV
2209 DO I =1,irun
2210 WU(I,L) = WU(I,L) + RHOKDZ(I,L) * ( U(I,L+1) - U(I,L) )
2211 WV(I,L) = WV(I,L) + RHOKDZ(I,L) * ( V(I,L+1) - V(I,L) )
2212 ENDDO
2213 ENDDO
2214 DO L =1,NLEVM1
2215 DO I =1,irun
2216 IF ( QQ(I,L) .GT. QQMIN ) THEN
2217 IBITSTB(I,L) = 1
2218 ELSE
2219 IBITSTB(I,L) = 0
2220 ENDIF
2221 IF( IBITSTB(I,L).EQ.1 ) FREQDG(I,L) = FREQDG(I,L) + aitr
2222 ENDDO
2223 ENDDO
1662f365b2 Andr*2224 do i = 1,irun
2225 qqcolmin(i) = qq(i,nlev)*0.1
2226 qqcolmax(i) = qq(i,nlev)
2227 levpbl(i) = nlev
2228 enddo
2229 DO L = nlev-1,1,-1
2230 DO I = 1,irun
2231 IF ( (qq(i,l).gt.qqcolmax(I)).and.(levpbl(i).eq.nlev))then
2232 qqcolmax(i) = qq(i,l)
2233 qqcolmin(i) = 0.1*qqcolmax(I)
2234 endif
2235 if((qq(i,l).lt.qqcolmin(i)).and.(levpbl(i).eq.nlev))
2236 1 levpbl(i)=l
2237 enddo
2238 enddo
2239 do i = 1,irun
2240 lp = levpbl(i)
2241 if(lp.lt.nlev)then
2242 pbldpth(I) = pbldpth(I) + ( (PLE(I,nlev+1)-PLE(I,Lp+2)) +
2243 1 ( (ple(i,lp+2)-ple(i,lp+1))*(qq(i,lp+1)-qqcolmin(i))
2244 2 / (qq(i,lp+1)-qq(i,lp)) ) ) * aitr
2245 else
2246 pbldpth(I) = pbldpth(I) + ( (PLE(I,nlev+1)-PLE(I,2)) +
2247 1 ( (ple(i,2)-ple(i,1))*(qq(i,1)-qqcolmin(i))
2248 2 / qq(i,1) ) ) * aitr
2249 endif
2250 enddo
2251 do i=1,irun
2252 sustar(i) = sustar(i) + aitr*ustar(i)
2253 enddo
2254 do i=1,irun
2255 sz0(i) = sz0(i) + aitr*z0(i)
2256 enddo
5317312052 Jean*2257 DO L =1,NLEV
2258 DO I =1,irun
2259 EU(I,L) = EU(I,L) + AITR*KM(I,L)
2260 ET(I,L) = ET(I,L) + AITR*KH(I,L)
2261 ENDDO
2262 ENDDO
1662f365b2 Andr*2263 DO I =1,irun
2264 scu(I) = scu(I) + AITR*cu(I)
2265 enddo
2266 DO I =1,irun
2267 sct(I) = sct(I) + AITR*ct(I)
2268 enddo
2269 IF(tprof) then
5317312052 Jean*2270 DO L =1,NLEVM1
2271 DO I =1,irun
2272 XLDIAG(I,L) = XLDIAG(I,L) + AITR*XL(I,L)
2273 ENDDO
2274 ENDDO
1662f365b2 Andr*2275 endif
2276 FIRST = .FALSE.
13fa5cb63c Jean*2277
1662f365b2 Andr*2278
13fa5cb63c Jean*2279
1662f365b2 Andr*2280 IF(ITER.EQ.ITRTRB)THEN
5317312052 Jean*2281 DO L =1,NLEVM1
2282 DO I =1,irun
2283 XLSAVE(I,L) = XL(I,L)
2284 ENDDO
1662f365b2 Andr*2285 ENDDO
2286 DO I = 1,irun
2287 CTSAVE(I) = CT(I)
2288 XXSAVE(I) = XX(I)
2289 YYSAVE(I) = YY(I)
2290 ZETASAVE(I) = ZETA(I)
2291 ENDDO
2292 ENDIF
2293
5317312052 Jean*2294 DO L =1,NLEV
2295 DO I =1,irun
2296 turbfcc(I,L) = turbfcc(I,L) + fcc(I,L) * aitr
2297 qliq(I,L) = qliq(I,L) + shl(I,L) * aitr
2298 ENDDO
2299 ENDDO
13fa5cb63c Jean*2300
1662f365b2 Andr*2301
13fa5cb63c Jean*2302
1662f365b2 Andr*2303 2000 CONTINUE
5317312052 Jean*2304 DO L =1,NLEV
2305 DO I =1,irun
2306 WU(I,L) = WU(I,L) * AITR
2307 WV(I,L) = WV(I,L) * AITR
2308 ENDDO
2309 ENDDO
13fa5cb63c Jean*2310
1662f365b2 Andr*2311 RETURN
2312 END
13fa5cb63c Jean*2313
2314
2315
1662f365b2 Andr*2316 SUBROUTINE SFCFLX(NN,VUS,VVS,VTHV1,VTHV2,VTH1,VTH2,VSH1,
2317 1 VSH2,VPK,VPKE,VPE,VZ0,IVWATER,VHS,
2318 2 VAHS,FIRST,LAST,N,IRUN,aitr,VRHO,VRHOZPK,VKH,VKM,
2319 3 VUSTAR,VXX,VYY,VCU,VCT,VRIB,VZETA,VWS,
2320 4 stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
2321 5 cp,rgas,undef,
2322 6 lwater, ivbitrib,
2323 7 VHZ,VPSIM,VAPSIM,VPSIG,VPSIHG,VTEMP,VDZETA,VDZ0,VDPSIM,
2324 8 VDPSIH,VZH,VXX0,VYY0,VAPSIHG,VRIB1,VWS1,VPSIH,VZETAL,
2325 9 VZ0L,VPSIH2,VH0,
2326 3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2327 4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
2328 5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
2329
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
2362
2363
2364
453e05dab3 Andr*2365 implicit none
1662f365b2 Andr*2366
453e05dab3 Andr*2367
2368 integer nn,n,irun
a456aa407c Andr*2369 _RL aitr,cp,rgas,undef
2370 _RL VUS(IRUN),VVS(IRUN),VTHV1(IRUN),VTHV2(IRUN)
2371 _RL VTH1(IRUN),VTH2(IRUN),VSH1(IRUN),VSH2(IRUN)
2372 _RL VPK(IRUN),VPKE(IRUN),VPE(IRUN)
2373 _RL VZ0(IRUN),VHS(IRUN),VAHS(IRUN)
453e05dab3 Andr*2374 integer IVWATER(IRUN)
2375 LOGICAL FIRST,LAST
a456aa407c Andr*2376 _RL VRHO(IRUN),VRHOZPK(IRUN)
2377 _RL VKM(IRUN),VKH(IRUN),VUSTAR(IRUN),VXX(IRUN)
2378 _RL VYY(IRUN),VCU(IRUN),VCT(IRUN),VRIB(IRUN)
2379 _RL VZETA(IRUN),VWS(IRUN)
2380 _RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
2381 _RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
453e05dab3 Andr*2382 LOGICAL LWATER
2383 integer IVBITRIB(irun)
a456aa407c Andr*2384 _RL VHZ(irun),VPSIM(irun),VAPSIM(irun),VPSIG(irun),VPSIHG(irun)
2385 _RL VTEMP(irun),VDZETA(irun),VDZ0(irun),VDPSIM(irun)
2386 _RL VDPSIH(irun),VZH(irun),VXX0(irun),VYY0(irun)
2387 _RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
2388 _RL VPSIH(irun),VZETAL(irun),VZ0L(irun),VPSIH2(irun),VH0(irun)
2389 _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
2390 _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
2391 _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
2392 _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
2393 _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
2394 _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
2395 _RL VDPSIMC(irun),VDPSIHC(irun)
453e05dab3 Andr*2396
2397
a456aa407c Andr*2398 _RL USTMX3,USTZ0S,Z0MIN,H0BYZ0,USTH0S,H0VEG,Z0VEGM,PRFAC
2399 _RL XPFAC,DIFSQT
1662f365b2 Andr*2400 PARAMETER ( USTMX3 = 0.0632456)
2401 PARAMETER ( USTZ0S = 0.2030325E-5)
2402 PARAMETER ( Z0MIN = USTZ0S/USTMX3)
2403 PARAMETER ( H0BYZ0 = 30.0 )
2404 PARAMETER ( USTH0S = H0BYZ0*USTZ0S )
5317312052 Jean*2405 PARAMETER ( H0VEG = 0.01 )
1662f365b2 Andr*2406 PARAMETER ( Z0VEGM = 0.005 )
2407 PARAMETER ( PRFAC = 0.595864 )
5317312052 Jean*2408 PARAMETER ( XPFAC = .55 )
1662f365b2 Andr*2409 PARAMETER ( DIFSQT = 3.872983E-3)
2410
a456aa407c Andr*2411 _RL psihdiag(irun),psimdiag(irun)
2412 _RL getcon,vk,rvk,vk2,bmdl
453e05dab3 Andr*2413 integer iwater,itype
2414 integer i,iter
13fa5cb63c Jean*2415
1662f365b2 Andr*2416 vk = getcon('VON KARMAN')
2417 rvk = 1./vk
2418 vk2 = vk*vk
2419 BMDL = VK * XPFAC * PRFAC / DIFSQT
2420
2421
13fa5cb63c Jean*2422
1662f365b2 Andr*2423 DO 9000 I = 1,IRUN
2424 VWS(I) = VUS(I) * VUS(I) + VVS(I) * VVS(I)
2425 IF ( VWS(I) .LE. 1.e-4) VWS(I) = 1.e-4
2426 VRIB(I) = ( CP * (VPKE(I)-VPK(I)) ) *
2427 1 (VTHV1(I) - VTHV2(I)) / VWS(I)
2428 VWS(I) = SQRT( VWS(I) )
2429 9000 CONTINUE
13fa5cb63c Jean*2430
1662f365b2 Andr*2431
2432
13fa5cb63c Jean*2433
1662f365b2 Andr*2434 IF (.NOT. FIRST) GO TO 100
13fa5cb63c Jean*2435
1662f365b2 Andr*2436 IWATER = 0
2437 DO 9002 I = 1,IRUN
2438 IF (IVWATER(I).EQ.1) IWATER = IWATER + 1
2439 9002 CONTINUE
2440 LWATER = .FALSE.
2441 IF(IWATER.GE.1)LWATER = .TRUE.
13fa5cb63c Jean*2442
1662f365b2 Andr*2443 IF(LWATER)THEN
2444 DO 9004 I = 1,IRUN
2445 IF (IVWATER(I).EQ.1) VZ0(I) = 0.0003
2446 9004 CONTINUE
2447 ENDIF
2448 do i = 1,irun
2449 vh0(i) = h0byz0 * vz0(i)
2450 if(vz0(i).ge.z0vegm)vh0(i) = h0veg
2451 enddo
2452
2453
13fa5cb63c Jean*2454
1662f365b2 Andr*2455 DO 9006 I = 1,IRUN
2456 VHZ(I) = VHS(I) / VZ0(I)
2457 VPSIM(I) = LOG( VHZ(I) )
2458 VAPSIM(I) = 1. / VPSIM(I)
2459 VCU(I) = VK * VAPSIM(I)
2460 VUSTAR(I) = VCU(I) * VWS(I)
13fa5cb63c Jean*2461
1662f365b2 Andr*2462 VPSIG(I) = VH0(I) * VUSTAR(I) - USTH0S
2463 if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2464 VPSIG(I) = SQRT( VPSIG(I) )
2465 VPSIG(I) = BMDL * VPSIG(I)
2466 VPSIHG(I) = VPSIM(I) + VPSIG(I)
2467 9006 CONTINUE
13fa5cb63c Jean*2468
1662f365b2 Andr*2469
13fa5cb63c Jean*2470
1662f365b2 Andr*2471 IF(LWATER)THEN
2472 DO 9008 I = 1,IRUN
2473 VTEMP(I) = 0.
2474 9008 CONTINUE
2475 CALL LINADJ(NN,VRIB,VRIB,VWS,
2476 1 VWS,VZ0,VUSTAR,IVWATER,
2477 2 VAPSIM,VTEMP,VTEMP,
2478 3 VTEMP,VTEMP,VTEMP,
2479 4 VTEMP,VTEMP,1,.TRUE.,IRUN,VDZETA,
2480 5 VDZ0,VDPSIM,VDPSIH,
2481 6 IVBITRIB,
2482 3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2483 4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
2484 5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
2485 DO 9010 I = 1,IRUN
2486 IF ( IVWATER(I).EQ.1 ) THEN
2487 VCU(I) = VCU(I) * (1. - VDPSIM(I)*VAPSIM(I))
2488 VZ0(I) = VZ0(I) + VDZ0(I)
2489 ENDIF
2490 IF ( IVWATER(I).EQ.1) THEN
2491 IF ( VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
2492 vh0(i) = h0byz0 * vz0(i)
2493 VPSIG(I) = VH0(I) * VCU(I) * VWS(I) - USTH0S
2494 if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2495 VPSIG(I) = SQRT( VPSIG(I) )
2496 VPSIG(I) = BMDL * VPSIG(I)
2497 VPSIHG(I) = VPSIM(I) + VDPSIH(I) + VPSIG(I)
2498 ENDIF
2499 9010 CONTINUE
2500 ENDIF
13fa5cb63c Jean*2501
1662f365b2 Andr*2502
13fa5cb63c Jean*2503
1662f365b2 Andr*2504 DO 9012 I = 1,IRUN
2505 VZETA(I) = VK2 * VRIB(I) / (VCU(I) * VCU(I) * VPSIHG(I))
2506 9012 CONTINUE
13fa5cb63c Jean*2507
1662f365b2 Andr*2508
13fa5cb63c Jean*2509
1662f365b2 Andr*2510 DO 9014 I = 1,IRUN
2511 VZH(I) = VZ0(I) * VAHS(I)
2512 9014 CONTINUE
2513 CALL PSI (VZETA,VZH,VPSIM,
2514 1 VTEMP,IRUN,VXX,VXX0,VYY,
2515 2 VYY0,2)
2516 DO 9016 I = 1,IRUN
2517 VCU(I) = VK / VPSIM(I)
2518 VPSIG(I) = VH0(I) * VCU(I) * VWS(I) - USTH0S
2519 if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2520 VPSIG(I) = SQRT(VPSIG(I))
2521 VPSIG(I) = BMDL * VPSIG(I)
2522 VPSIHG(I) = VPSIM(I) + VPSIG(I)
2523 VZETA(I) = VK2 * VRIB(I) / (VCU(I) * VCU(I) * VPSIHG(I))
2524 9016 CONTINUE
13fa5cb63c Jean*2525
1662f365b2 Andr*2526 IF(LWATER)THEN
2527
2528 DO 9018 I = 1,IRUN
2529 IF (IVWATER(I).EQ.1) VUSTAR(I) = VCU(I) * VWS(I)
2530 9018 CONTINUE
2531 CALL ZCSUB ( VUSTAR,VHZ,IVWATER,.FALSE.,IRUN,VTEMP)
2532 DO 9020 I = 1,IRUN
2533 IF (IVWATER(I).EQ.1 ) then
2534 VZ0(I) = VTEMP(I)
2535 IF ( VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
2536 vh0(i) = h0byz0 * vz0(i)
2537 endif
2538 9020 CONTINUE
2539 ENDIF
13fa5cb63c Jean*2540
1662f365b2 Andr*2541 GO TO 125
13fa5cb63c Jean*2542
1662f365b2 Andr*2543
13fa5cb63c Jean*2544
1662f365b2 Andr*2545 100 CONTINUE
2546
2547 CALL LINADJ(NN,VRIB1,VRIB,VWS1,
2548 1 VWS,VZ0,VUSTAR,IVWATER,
2549 2 VAPSIM,VAPSIHG,VPSIH,
2550 3 VPSIG,VXX,VXX0,
2551 4 VYY,VYY0,2,LWATER,IRUN,VDZETA,
2552 5 VDZ0,VDPSIM,VDPSIH,
2553 6 IVBITRIB,
2554 3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2555 4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
2556 5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
13fa5cb63c Jean*2557
1662f365b2 Andr*2558 DO 9022 I = 1,IRUN
2559 VZETA(I) = VZETA(I) + VZETAL(I) * VDZETA(I)
2560 IF (IVBITRIB(I).EQ.1 )VZETA(I) =
2561 1 VPSIM(I) * VPSIM(I) * VRIB(I) * VCT(I) * RVK
2562 9022 CONTINUE
13fa5cb63c Jean*2563
1662f365b2 Andr*2564 IF ( LWATER ) THEN
2565 DO 9024 I = 1,IRUN
2566 IF (IVWATER(I).EQ.1 ) then
2567 VZ0(I) = VZ0(I) + VZ0L(I) * VDZ0(I)
2568 IF (VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
2569 vh0(i) = h0byz0 * vz0(i)
2570 endif
2571 9024 CONTINUE
2572 ENDIF
13fa5cb63c Jean*2573
1662f365b2 Andr*2574 125 CONTINUE
13fa5cb63c Jean*2575
1662f365b2 Andr*2576
2577
13fa5cb63c Jean*2578
1662f365b2 Andr*2579 DO 200 ITER = 1,N
2580 DO 9026 I = 1,IRUN
2581 VZH(I) = VZ0(I) * VAHS(I)
2582 9026 CONTINUE
2583 CALL PSI (VZETA,VZH,VPSIM,
2584 1 VPSIH,IRUN,VXX,VXX0,VYY,
2585 2 VYY0,1)
2586 DO 9028 I = 1,IRUN
2587 VCU(I) = VK / VPSIM(I)
2588 VUSTAR(I) = VCU(I) * VWS(I)
13fa5cb63c Jean*2589
1662f365b2 Andr*2590 VPSIG(I) = VH0(I) * VUSTAR(I) - USTH0S
2591 if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2592 VPSIG(I) = SQRT(VPSIG(I))
2593 VPSIG(I) = BMDL * VPSIG(I)
2594 VPSIHG(I) = VPSIH(I) + VPSIG(I)
13fa5cb63c Jean*2595
1662f365b2 Andr*2596
13fa5cb63c Jean*2597
1662f365b2 Andr*2598 VAPSIM(I) = VCU(I) * RVK
2599 VAPSIHG(I) = 1. / VPSIHG(I)
2600 VRIB1(I) = VAPSIM(I) * VAPSIM(I) * VPSIHG(I) * VZETA(I)
2601 9028 CONTINUE
13fa5cb63c Jean*2602
1662f365b2 Andr*2603 ITYPE = 3
2604 IF(ITER.EQ.N) ITYPE = 4
2605 IF( (ITYPE.EQ.4) .AND. (.NOT.LAST) ) ITYPE = 5
13fa5cb63c Jean*2606
1662f365b2 Andr*2607 CALL LINADJ(NN,VRIB1,VRIB,VWS,
2608 1 VWS,VZ0,VUSTAR,IVWATER,
2609 2 VAPSIM,VAPSIHG,VPSIH,
2610 3 VPSIG,VXX,VXX0,
2611 4 VYY,VYY0,ITYPE,LWATER,IRUN,VDZETA,
2612 5 VDZ0,VDPSIM,VDPSIH,
2613 6 IVBITRIB,
2614 3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2615 4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
2616 5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
13fa5cb63c Jean*2617
1662f365b2 Andr*2618
13fa5cb63c Jean*2619
1662f365b2 Andr*2620 IF (ITYPE.EQ.5) THEN
2621 DO 9030 I = 1,IRUN
2622 VZETAL(I) = VZETA(I)
2623 VZ0L(I) = VZ0(I)
2624 9030 CONTINUE
2625 ENDIF
13fa5cb63c Jean*2626
1662f365b2 Andr*2627 DO 9032 I = 1,IRUN
2628 VZETA(I) = VZETA(I) * ( 1. + VDZETA(I) )
2629 IF (IVBITRIB(I).EQ.1 ) VZETA(I) =
2630 1 VPSIM(I) * VPSIM(I) * VRIB(I) * VAPSIHG(I)
2631 9032 CONTINUE
13fa5cb63c Jean*2632
1662f365b2 Andr*2633 IF ( LWATER ) THEN
2634 DO 9034 I = 1,IRUN
2635 IF (IVWATER(I).EQ.1 ) then
2636 VZ0(I) = VZ0(I) * ( 1. + VDZ0(I) )
2637 IF (VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
2638 vh0(i) = h0byz0 * vz0(i)
2639 endif
2640 9034 CONTINUE
2641 ENDIF
13fa5cb63c Jean*2642
1662f365b2 Andr*2643 IF ( ITER .EQ. N ) THEN
2644 DO 9036 I = 1,IRUN
2645 VPSIM(I) = VPSIM(I) + VDPSIM(I)
2646 VCU(I) = VK / VPSIM(I)
2647 VUSTAR(I) = VCU(I) * VWS(I)
13fa5cb63c Jean*2648
1662f365b2 Andr*2649 VPSIG(I) = VH0(I) * VUSTAR(I) - USTH0S
2650 if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2651 VPSIG(I) = SQRT(VPSIG(I))
2652 VPSIG(I) = BMDL * VPSIG(I)
2653 VPSIHG(I) = VPSIH(I) + VDPSIH(I) + VPSIG(I)
2654 VCT(I) = VK / VPSIHG(I)
2655 9036 CONTINUE
2656 ENDIF
13fa5cb63c Jean*2657
1662f365b2 Andr*2658
13fa5cb63c Jean*2659
1662f365b2 Andr*2660 IF (ITYPE.EQ.5) THEN
2661 DO 9038 I = 1,IRUN
2662 VRIB1(I) = VRIB(I)
2663 VWS1(I) = VWS(I)
2664 9038 CONTINUE
2665 ENDIF
13fa5cb63c Jean*2666
1662f365b2 Andr*2667 200 CONTINUE
13fa5cb63c Jean*2668
1662f365b2 Andr*2669
13fa5cb63c Jean*2670
1662f365b2 Andr*2671 IF (FIRST) THEN
2672 DO I = 1,IRUN
2673 VTEMP(I) = 10. * VAHS(I) * VZETA(I)
2674 VZH(I) = VZ0(I) * 0.1
2675 ENDDO
2676 CALL PSI (VTEMP,VZH,VHZ,
2677 1 VPSIH2,IRUN,VHZ,VHZ,VHZ,
2678 2 VHZ,3)
2679 DO I = 1,IRUN
2680 VTEMP(I) = ( VPSIH2(I) + VPSIG(I) ) / VPSIHG(I)
2681 VRHO(I) = VPKE(I)*( VTH2(I) + VTEMP(I) * (VTH1(I)-VTH2(I)) )
2682 VRHO(I) = VPE(I)*100. / ( RGAS * VRHO(I) )
2683 ENDDO
2684 ENDIF
13fa5cb63c Jean*2685
1662f365b2 Andr*2686
2687
2688
2689
13fa5cb63c Jean*2690
1662f365b2 Andr*2691 do i = 1,irun
2692 vtemp(i) = 2. * vahs(i) * vzeta(i)
2693 vzh(i) = vz0(i) * 0.5
2694 if(vz0(i).ge.2.)vzh(i) = 0.9
2695 enddo
2696 call psi(vtemp,vzh,psimdiag,psihdiag,irun,vhz,vhz,vhz,vhz,1)
2697 do i = 1,irun
2698 stu2m(i) = (psimdiag(i)/vpsim(i) * vus(i))
2699 stv2m(i) = (psimdiag(i)/vpsim(i) * vvs(i))
2700 stt2m(i) = ( (vth2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
2701 1 (vth1(i)-vth2(i))) ) * vpke(i)
2702 stq2m(i) = (vsh2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
2703 1 (vsh1(i)-vsh2(i)))
2704 if(vz0(i).ge.2.)then
2705 stu2m(i) = UNDEF
2706 stv2m(i) = UNDEF
2707 stt2m(i) = UNDEF
2708 stq2m(i) = UNDEF
2709 endif
2710 enddo
2711 do i = 1,irun
2712 vtemp(i) = 10. * vahs(i) * vzeta(i)
2713 vzh(i) = vz0(i) * 0.1
2714 enddo
2715 call psi(vtemp,vzh,psimdiag,psihdiag,irun,vhz,vhz,vhz,vhz,1)
2716 do i = 1,irun
2717 stu10m(i) = (psimdiag(i)/vpsim(i) * vus(i))
2718 stv10m(i) = (psimdiag(i)/vpsim(i) * vvs(i))
2719 stt10m(i) = ( (vth2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
2720 1 (vth1(i)-vth2(i))) ) * vpke(i)
2721 stq10m(i) = (vsh2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
2722 1 (vsh1(i)-vsh2(i)))
2723 enddo
13fa5cb63c Jean*2724
1662f365b2 Andr*2725
13fa5cb63c Jean*2726
1662f365b2 Andr*2727 DO 9044 I = 1,IRUN
2728 VRHOZPK(I) = VRHO(I) * VPKE(I)
2729 VKH(I) = VUSTAR(I) * VCT(I)
2730 VKM(I) = VUSTAR(I) * VCU(I)
2731 9044 CONTINUE
13fa5cb63c Jean*2732
1662f365b2 Andr*2733 RETURN
2734 END
13fa5cb63c Jean*2735
2736
2737
1662f365b2 Andr*2738 SUBROUTINE PHI(Z,PHIM,PHIH,IFLAG,N)
2739
2740
2741
2742
2743
2744
2745
2746
2747
2748
2749
2750
453e05dab3 Andr*2751 implicit none
2752
2753
2754 integer n,iflag
a456aa407c Andr*2755 _RL PHIM(N),PHIH(N),Z(N)
453e05dab3 Andr*2756
2757
2758 integer I1(N),I2(N)
a456aa407c Andr*2759 _RL ZSTAR(N),E1(N),E2(N),TEMP1(N)
13fa5cb63c Jean*2760
a456aa407c Andr*2761 _RL PHIM0(385),ZLINM1(75),ZLINM2(75),ZLINM3(36)
2762 _RL ZLOGM1(74),ZLOGM2(75),ZLOGM3(50)
2763 _RL PHIH0(385),ZLINH1(75),ZLINH2(75),ZLINH3(36)
2764 _RL ZLOGH1(74),ZLOGH2(75),ZLOGH3(50)
1662f365b2 Andr*2765 EQUIVALENCE (PHIM0(1),ZLINM1(1)),(PHIM0(76),ZLINM2(1))
2766 EQUIVALENCE (PHIM0(151),ZLINM3(1))
2767 EQUIVALENCE (PHIM0(187),ZLOGM1(1)),(PHIM0(261),ZLOGM2(1))
2768 EQUIVALENCE (PHIM0(336),ZLOGM3(1))
2769 EQUIVALENCE (PHIH0(1),ZLINH1(1)),(PHIH0(76),ZLINH2(1))
2770 EQUIVALENCE (PHIH0(151),ZLINH3(1))
2771 EQUIVALENCE (PHIH0(187),ZLOGH1(1)),(PHIH0(261),ZLOGH2(1))
2772 EQUIVALENCE (PHIH0(336),ZLOGH3(1))
13fa5cb63c Jean*2773
1662f365b2 Andr*2774 DATA ZLOGM1/
13fa5cb63c Jean*2775 & 0.697894,0.678839,0.659598,0.640260,
2776 & 0.620910,0.601628,0.582486,0.563550,0.544877,
2777 & 0.526519,0.508516,0.490903,0.473708,0.456951,
2778 & 0.440649,0.424812,0.409446,0.394553,0.380133,
2779 & 0.366182,0.352695,0.339664,0.327082,0.314938,
2780 & 0.303222,0.291923,0.281029,0.270528,0.260409,
2781 & 0.250659,0.241267,0.232221,0.223509,0.215119,
2782 & 0.207041,0.199264,0.191776,0.184568,0.177628,
2783 & 0.170949,0.164519,0.158331,0.152374,0.146641,
2784 & 0.141123,0.135813,0.130702,0.125783,0.121048,
2785 & 0.116492,0.112107,0.107887,0.103826,0.0999177,
2786 & 0.0961563,0.0925364,0.0890528,0.0857003,0.0824739,
2787 & 0.0793690,0.0763810,0.0735054,0.0707380,0.0680749,
2788 & 0.0655120,0.0630455,0.0606720,0.0583877,0.0561895,
2789 & 0.0540740,0.0520382,0.0500790,0.0481936,0.0463791/
1662f365b2 Andr*2790 DATA ZLOGM2/
13fa5cb63c Jean*2791 & 0.0446330,0.0429526,0.0413355,0.0397792,0.0382816,
2792 & 0.0368403,0.0354533,0.0341185,0.0328340,0.0315978,
2793 & 0.0304081,0.0292633,0.0281616,0.0271013,0.0260809,
2794 & 0.0250990,0.0241540,0.0232447,0.0223695,0.0215273,
2795 & 0.0207168,0.0199369,0.0191862,0.0184639,0.0177687,
2796 & 0.0170998,0.0164560,0.0158364,0.0152402,0.0146664,
2797 & 0.0141142,0.0135828,0.0130714,0.0125793,0.0121057,
2798 & 0.0116499,0.0112113,0.0107892,0.0103830,0.999210E-2,
2799 & 0.961590E-2,0.925387E-2,0.890547E-2,0.857018E-2,0.824752E-2,
2800 & 0.793701E-2,0.763818E-2,0.735061E-2,0.707386E-2,0.680754E-2,
2801 & 0.655124E-2,0.630459E-2,0.606722E-2,0.583880E-2,0.561897E-2,
2802 & 0.540742E-2,0.520383E-2,0.500791E-2,0.481937E-2,0.463792E-2,
2803 & 0.446331E-2,0.429527E-2,0.413355E-2,0.397793E-2,0.382816E-2,
2804 & 0.368403E-2,0.354533E-2,0.341185E-2,0.328340E-2,0.315978E-2,
2805 & 0.304082E-2,0.292633E-2,0.281616E-2,0.271013E-2,0.260809E-2/
1662f365b2 Andr*2806 DATA ZLOGM3/
13fa5cb63c Jean*2807 & 0.250990E-2,0.241541E-2,0.232447E-2,0.223695E-2,0.215273E-2,
2808 & 0.207168E-2,0.199369E-2,0.191862E-2,0.184639E-2,0.177687E-2,
2809 & 0.170998E-2,0.164560E-2,0.158364E-2,0.152402E-2,0.146664E-2,
2810 & 0.141142E-2,0.135828E-2,0.130714E-2,0.125793E-2,0.121057E-2,
2811 & 0.116499E-2,0.112113E-2,0.107892E-2,0.103830E-2,0.999210E-3,
2812 & 0.961590E-3,0.925387E-3,0.890547E-3,0.857018E-3,0.824752E-3,
2813 & 0.793701E-3,0.763818E-3,0.735061E-3,0.707386E-3,0.680754E-3,
2814 & 0.655124E-3,0.630459E-3,0.606722E-3,0.583880E-3,0.561897E-3,
2815 & 0.540742E-3,0.520383E-3,0.500791E-3,0.481937E-3,0.463792E-3,
2816 & 0.446331E-3,0.429527E-3,0.413355E-3,0.397793E-3,0.382816E-3/
1662f365b2 Andr*2817 DATA ZLOGH1/
13fa5cb63c Jean*2818 & 0.640529,0.623728,0.606937,0.590199,
2819 & 0.573552,0.557032,0.540672,0.524504,0.508553,
2820 & 0.492843,0.477397,0.462232,0.447365,0.432809,
2821 & 0.418574,0.404670,0.391103,0.377878,0.364999,
2822 & 0.352468,0.340284,0.328447,0.316954,0.305804,
2823 & 0.294992,0.284514,0.274364,0.264538,0.255028,
2824 & 0.245829,0.236933,0.228335,0.220026,0.211999,
2825 & 0.204247,0.196762,0.189537,0.182564,0.175837,
2826 & 0.169347,0.163088,0.157051,0.151231,0.145620,
2827 & 0.140211,0.134998,0.129974,0.125133,0.120469,
2828 & 0.115975,0.111645,0.107475,0.103458,0.995895E-1,
2829 & 0.958635E-1,0.922753E-1,0.888199E-1,0.854925E-1,0.822886E-1,
2830 & 0.792037E-1,0.762336E-1,0.733739E-1,0.706208E-1,0.679704E-1,
2831 & 0.654188E-1,0.629625E-1,0.605979E-1,0.583217E-1,0.561306E-1,
2832 & 0.540215E-1,0.519914E-1,0.500373E-1,0.481564E-1,0.463460E-1/
1662f365b2 Andr*2833 DATA ZLOGH2/
13fa5cb63c Jean*2834 & 0.446034E-1,0.429263E-1,0.413120E-1,0.397583E-1,0.382629E-1,
2835 & 0.368237E-1,0.354385E-1,0.341053E-1,0.328222E-1,0.315873E-1,
2836 & 0.303988E-1,0.292550E-1,0.281541E-1,0.270947E-1,0.260750E-1,
2837 & 0.250937E-1,0.241494E-1,0.232405E-1,0.223658E-1,0.215240E-1,
2838 & 0.207139E-1,0.199342E-1,0.191839E-1,0.184618E-1,0.177669E-1,
2839 & 0.170981E-1,0.164545E-1,0.158351E-1,0.152390E-1,0.146653E-1,
2840 & 0.141133E-1,0.135820E-1,0.130707E-1,0.125786E-1,0.121051E-1,
2841 & 0.116494E-1,0.112108E-1,0.107888E-1,0.103826E-1,0.999177E-2,
2842 & 0.961561E-2,0.925360E-2,0.890523E-2,0.856997E-2,0.824733E-2,
2843 & 0.793684E-2,0.763803E-2,0.735048E-2,0.707375E-2,0.680743E-2,
2844 & 0.655114E-2,0.630450E-2,0.606715E-2,0.583873E-2,0.561891E-2,
2845 & 0.540737E-2,0.520379E-2,0.500787E-2,0.481933E-2,0.463789E-2,
2846 & 0.446328E-2,0.429524E-2,0.413353E-2,0.397790E-2,0.382814E-2,
2847 & 0.368401E-2,0.354532E-2,0.341184E-2,0.328338E-2,0.315977E-2,
2848 & 0.304081E-2,0.292632E-2,0.281615E-2,0.271012E-2,0.260809E-2/
1662f365b2 Andr*2849 DATA ZLOGH3/
13fa5cb63c Jean*2850 & 0.250990E-2,0.241540E-2,0.232446E-2,0.223695E-2,0.215273E-2,
2851 & 0.207168E-2,0.199368E-2,0.191862E-2,0.184639E-2,0.177687E-2,
2852 & 0.170997E-2,0.164559E-2,0.158364E-2,0.152402E-2,0.146664E-2,
2853 & 0.141142E-2,0.135828E-2,0.130714E-2,0.125793E-2,0.121057E-2,
2854 & 0.116499E-2,0.112113E-2,0.107892E-2,0.103830E-2,0.999209E-3,
2855 & 0.961590E-3,0.925387E-3,0.890546E-3,0.857018E-3,0.824752E-3,
2856 & 0.793700E-3,0.763818E-3,0.735061E-3,0.707386E-3,0.680754E-3,
2857 & 0.655124E-3,0.630459E-3,0.606722E-3,0.583880E-3,0.561897E-3,
2858 & 0.540742E-3,0.520383E-3,0.500791E-3,0.481937E-3,0.463792E-3,
2859 & 0.446331E-3,0.429527E-3,0.413355E-3,0.397793E-3,0.382816E-3/
2860
1662f365b2 Andr*2861 DATA ZLINM1/
2862 & 0.964508,0.962277,0.960062,0.957863,0.955680,
2863 & 0.953512,0.951359,0.949222,0.947100,0.944992,
2864 & 0.942899,0.940821,0.938758,0.936709,0.934673,
2865 & 0.932652,0.930645,0.928652,0.926672,0.924706,
2866 & 0.922753,0.920813,0.918886,0.916973,0.915072,
2867 & 0.913184,0.911308,0.909445,0.907594,0.905756,
2868 & 0.903930,0.902115,0.900313,0.898522,0.896743,
2869 & 0.894975,0.893219,0.891475,0.889741,0.888019,
2870 & 0.886307,0.884607,0.882917,0.881238,0.879569,
2871 & 0.877911,0.876264,0.874626,0.872999,0.871382,
2872 & 0.869775,0.868178,0.866591,0.865013,0.863445,
2873 & 0.861887,0.860338,0.858798,0.857268,0.855747,
2874 & 0.854235,0.852732,0.851238,0.849753,0.848277,
2875 & 0.846809,0.845350,0.843900,0.842458,0.841025,
2876 & 0.839599,0.838182,0.836774,0.835373,0.833980/
2877 DATA ZLINM2/
2878 & 0.832596,0.831219,0.829850,0.828489,0.827136,
2879 & 0.825790,0.824451,0.823121,0.821797,0.820481,
2880 & 0.819173,0.817871,0.816577,0.815289,0.814009,
2881 & 0.812736,0.811470,0.810210,0.808958,0.807712,
2882 & 0.806473,0.805240,0.804015,0.802795,0.801582,
2883 & 0.800376,0.799176,0.797982,0.796794,0.795613,
2884 & 0.794438,0.793269,0.792106,0.790949,0.789798,
2885 & 0.788652,0.787513,0.786380,0.785252,0.784130,
2886 & 0.783014,0.781903,0.780798,0.779698,0.778604,
2887 & 0.777516,0.776432,0.775354,0.774282,0.773215,
2888 & 0.772153,0.771096,0.770044,0.768998,0.767956,
2889 & 0.766920,0.765888,0.764862,0.763840,0.762824,
2890 & 0.761812,0.760805,0.759803,0.758805,0.757813,
2891 & 0.756824,0.755841,0.754862,0.753888,0.752918,
2892 & 0.751953,0.750992,0.750035,0.749083,0.748136/
2893 DATA ZLINM3/
2894 & 0.747192,0.746253,0.745318,0.744388,0.743462,
2895 & 0.742539,0.741621,0.740707,0.739798,0.738892,
2896 & 0.737990,0.737092,0.736198,0.735308,0.734423,
2897 & 0.733540,0.732662,0.731788,0.730917,0.730050,
2898 & 0.729187,0.728328,0.727472,0.726620,0.725772,
2899 & 0.724927,0.724086,0.723248,0.722414,0.721584,
2900 & 0.720757,0.719933,0.719113,0.718296,0.717483,
2901 & 0.716673/
2902 DATA ZLINH1/
2903 & 0.936397,0.932809,0.929287,0.925827,0.922429,
2904 & 0.919089,0.915806,0.912579,0.909405,0.906284,
2905 & 0.903212,0.900189,0.897214,0.894284,0.891399,
2906 & 0.888558,0.885759,0.883001,0.880283,0.877603,
2907 & 0.874962,0.872357,0.869788,0.867255,0.864755,
2908 & 0.862288,0.859854,0.857452,0.855081,0.852739,
2909 & 0.850427,0.848144,0.845889,0.843662,0.841461,
2910 & 0.839287,0.837138,0.835014,0.832915,0.830841,
2911 & 0.828789,0.826761,0.824755,0.822772,0.820810,
2912 & 0.818869,0.816949,0.815050,0.813170,0.811310,
2913 & 0.809470,0.807648,0.805845,0.804060,0.802293,
2914 & 0.800543,0.798811,0.797095,0.795396,0.793714,
2915 & 0.792047,0.790396,0.788761,0.787141,0.785535,
2916 & 0.783945,0.782369,0.780807,0.779259,0.777724,
2917 & 0.776204,0.774696,0.773202,0.771720,0.770251/
2918 DATA ZLINH2/
2919 & 0.768795,0.767351,0.765919,0.764499,0.763091,
2920 & 0.761694,0.760309,0.758935,0.757571,0.756219,
2921 & 0.754878,0.753547,0.752226,0.750916,0.749616,
2922 & 0.748326,0.747045,0.745775,0.744514,0.743262,
2923 & 0.742020,0.740787,0.739563,0.738348,0.737141,
2924 & 0.735944,0.734755,0.733574,0.732402,0.731238,
2925 & 0.730083,0.728935,0.727795,0.726664,0.725539,
2926 & 0.724423,0.723314,0.722213,0.721119,0.720032,
2927 & 0.718952,0.717880,0.716815,0.715756,0.714704,
2928 & 0.713660,0.712621,0.711590,0.710565,0.709547,
2929 & 0.708534,0.707529,0.706529,0.705536,0.704549,
2930 & 0.703567,0.702592,0.701623,0.700660,0.699702,
2931 & 0.698750,0.697804,0.696863,0.695928,0.694998,
2932 & 0.694074,0.693155,0.692241,0.691333,0.690430,
2933 & 0.689532,0.688639,0.687751,0.686868,0.685990/
2934 DATA ZLINH3/
2935 & 0.685117,0.684249,0.683386,0.682527,0.681673,
2936 & 0.680824,0.679979,0.679139,0.678303,0.677472,
2937 & 0.676645,0.675823,0.675005,0.674191,0.673381,
2938 & 0.672576,0.671775,0.670978,0.670185,0.669396,
2939 & 0.668611,0.667830,0.667054,0.666281,0.665512,
2940 & 0.664746,0.663985,0.663227,0.662473,0.661723,
2941 & 0.660977,0.660234,0.659495,0.658759,0.658027,
2942 & 0.657298/
453e05dab3 Andr*2943
2944 integer ibit1,ibit2,i
13fa5cb63c Jean*2945
1662f365b2 Andr*2946 IBIT1 = 0
2947 IBIT2 = 0
13fa5cb63c Jean*2948
1662f365b2 Andr*2949 DO 9000 I = 1,N
2950 IF(Z(I).GE.0.15)IBIT1 = IBIT1 + 1
2951 IF(Z(I).GT.2. )IBIT2 = IBIT2 + 1
2952 9000 CONTINUE
13fa5cb63c Jean*2953
1662f365b2 Andr*2954 IF( IBIT1 .LE. 0 ) GO TO 200
13fa5cb63c Jean*2955
1662f365b2 Andr*2956 DO 9002 I = 1,N
2957 ZSTAR(I) = 100. * Z(I) - 14.
2958 9002 CONTINUE
13fa5cb63c Jean*2959
1662f365b2 Andr*2960 IF( IBIT2 .LE. 0 ) GO TO 60
2961 DO 9004 I = 1,N
2962 TEMP1(I) = Z(I)*0.5
2963 IF( Z(I) .LE. 2. )TEMP1(I) = 1.
2964 TEMP1(I) = LOG10(TEMP1(I))
2965 TEMP1(I) = (TEMP1(I) + 9.3) * 20.
2966 IF( Z(I) .GT. 2. ) ZSTAR(I) = TEMP1(I)
2967 IF( Z(I).GT.1.78e10 ) ZSTAR(I) = 384.9999
2968 9004 CONTINUE
13fa5cb63c Jean*2969
1662f365b2 Andr*2970 60 CONTINUE
13fa5cb63c Jean*2971
1662f365b2 Andr*2972 DO 9006 I = 1,N
2973 I1(I) = ZSTAR(I)
2974 I2(I) = I1(I) + 1
2975 TEMP1(I) = ZSTAR(I) - I1(I)
13fa5cb63c Jean*2976
1662f365b2 Andr*2977 9006 CONTINUE
13fa5cb63c Jean*2978
1662f365b2 Andr*2979 IF( IFLAG .GT. 2 ) GO TO 100
2980 DO 9008 I = 1,N
2981 if( z(i).ge.0.15 ) then
2982 E1(I) = PHIM0( I1(I) )
2983 E2(I) = PHIM0( I2(I) )
2984 PHIM(I) = TEMP1(I) * ( E2(I)-E1(I) )
2985 PHIM(I) = PHIM(I) + E1(I)
2986 endif
2987 9008 CONTINUE
2988
2989 100 CONTINUE
13fa5cb63c Jean*2990
1662f365b2 Andr*2991 IF( IFLAG .EQ. 2 ) GO TO 200
2992 DO 9010 I = 1,N
2993 if( z(i).ge.0.15 ) then
2994 E1(I) = PHIH0( I1(I) )
2995 E2(I) = PHIH0( I2(I) )
2996 PHIH(I) = TEMP1(I) * ( E2(I)-E1(I) )
2997 PHIH(I) = PHIH(I) + E1(I)
2998 endif
2999 9010 CONTINUE
3000
3001 200 CONTINUE
3002 IF( IBIT1 .GE. N ) GO TO 500
13fa5cb63c Jean*3003
1662f365b2 Andr*3004 DO 9012 I = 1,N
3005 ZSTAR(I) = -Z(I)
3006 9012 CONTINUE
13fa5cb63c Jean*3007
1662f365b2 Andr*3008 IF( IFLAG .GT. 2 ) GO TO 300
3009 DO 9014 I = 1,N
3010 IF( Z(I) .LT. 0.15 ) PHIM(I) = 1. + ZSTAR(I)
3011 2 *(0.25+ZSTAR(I)*(0.09375+ZSTAR(I)*
3012 3 (0.03125+0.00732422 * ZSTAR(I))))
3013 9014 CONTINUE
13fa5cb63c Jean*3014
1662f365b2 Andr*3015 300 CONTINUE
3016 IF( IFLAG .EQ. 2 ) GO TO 500
3017 DO 9016 I = 1,N
3018 IF( Z(I) .LT. 0.15 ) THEN
3019 PHIH(I) =1.+ Z(I) * (0.5+ZSTAR(I)*(0.375+ZSTAR(I)*
3020 1 (0.5+ZSTAR(I)*(0.8203125+ZSTAR(I)*
3021 2 (1.5+2.93262*ZSTAR(I))))))
3022 PHIH(I) = 1. / PHIH(I)
3023 ENDIF
3024 9016 CONTINUE
13fa5cb63c Jean*3025
1662f365b2 Andr*3026 500 CONTINUE
3027 RETURN
3028 END
13fa5cb63c Jean*3029
3030
3031
1662f365b2 Andr*3032 SUBROUTINE PSI(VZZ,VZH,VPSIM,VPSIH,IRUN,VX,VXS,VY,VYS,IFLAG)
3033
3034
3035
3036
3037
3038
3039
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049
3050
453e05dab3 Andr*3051 implicit none
3052
3053
3054 integer irun,iflag
a456aa407c Andr*3055 _RL VZZ(IRUN),VZH(IRUN),VPSIM(IRUN),VPSIH(IRUN),
453e05dab3 Andr*3056 1 VX(IRUN),VXS(IRUN),VY(IRUN),VYS(IRUN)
5317312052 Jean*3057
453e05dab3 Andr*3058
a456aa407c Andr*3059 _RL ZWM,RZWM,Z0M,ZCM,RZCM,CM1,CM2,CM6,CM7,CM8ARG,YCM
1662f365b2 Andr*3060 PARAMETER ( ZWM = 1. )
3061 PARAMETER ( RZWM = 1./ZWM )
3062 PARAMETER ( Z0M = 0.2 )
3063 PARAMETER ( ZCM = 42. )
3064 PARAMETER ( RZCM = 1./ZCM )
3065 PARAMETER ( CM1 = 1./126. )
3066 PARAMETER ( CM2 = 1./(6.*CM1) )
3067 PARAMETER ( CM6 = 6. / ( 1. + 6.*CM1 ) )
3068 PARAMETER ( CM7 = CM2 + ZWM )
3069 PARAMETER ( CM8ARG = CM7*ZCM*RZWM / (CM2+ZCM) )
3070 PARAMETER ( YCM = 6. / ( 1. + 6.*CM1*ZCM ) )
3071
453e05dab3 Andr*3072 integer INTSTB(irun),INTZ0(irun)
a456aa407c Andr*3073 _RL ZZ0(irun),Z(irun),Z2(irun),Z1(irun),Z0(irun)
3074 _RL X0(irun),X1(irun),Y0(irun),Y1(irun)
3075 _RL PSI2(irun),TEMP(irun)
3076 _RL HZ(irun),ARG0(irun),ARG1(irun),DX(irun)
3077 _RL X0NUM(irun),X1NUM(irun),X0DEN(irun)
3078 _RL X1DEN(irun),Y1DEN(irun),Z2ZWM(irun)
3079 _RL cm3,cm4,cm5,cm8
471f7f004b Andr*3080 integer ibit,indx
453e05dab3 Andr*3081 integer i
13fa5cb63c Jean*3082
1662f365b2 Andr*3083 CM3 = sqrt( 0.2/CM1-0.01 )
3084 CM4 = 1./CM3
3085 CM5 = (10.-CM1) / (10.*CM1*CM3)
3086 CM8 = 6. * LOG(CM8ARG)
13fa5cb63c Jean*3087
1662f365b2 Andr*3088 DO 9000 I = 1,IRUN
3089 VPSIM(I) = 0.
3090 VPSIH(I) = 0.
3091 VX(I) = 0.
3092 VXS(I) = 0.
3093 VY(I) = 0.
3094 VYS(I) = 0.
3095 ZZ0(I) = VZH(I)*VZZ(I)
3096 9000 CONTINUE
3097 IBIT = 0
3098 DO 9122 I = 1,IRUN
3099 IF(VZZ(I).LE.-1.e-7)IBIT = IBIT + 1
3100 9122 CONTINUE
3101 DO 9022 I = 1,IRUN
3102 IF(VZZ(I).LE.-1.e-7)THEN
3103 INTSTB(I) = 1
3104 ELSE
3105 INTSTB(I) = 0
3106 ENDIF
3107 9022 CONTINUE
13fa5cb63c Jean*3108
1662f365b2 Andr*3109
3110
3111
13fa5cb63c Jean*3112
1662f365b2 Andr*3113 IF(IBIT.LE.0) GO TO 100
13fa5cb63c Jean*3114
471f7f004b Andr*3115 indx = 0
1662f365b2 Andr*3116 DO 9002 I = 1,IRUN
3117 IF (INTSTB(I).EQ.1)THEN
471f7f004b Andr*3118 indx = indx + 1
3119 Z(indx) = VZZ(I)
3120 Z0(indx) = ZZ0(I)
1662f365b2 Andr*3121 ENDIF
3122 9002 CONTINUE
13fa5cb63c Jean*3123
1662f365b2 Andr*3124 DO 9004 I = 1,IBIT
3125 Z(I) = -18. * Z(I)
3126 Z0(I) = -18. * Z0(I)
3127 9004 CONTINUE
5317312052 Jean*3128
1662f365b2 Andr*3129 CALL PHI( Z,X1,Y1,IFLAG,IBIT )
3130 CALL PHI( Z0,X0,Y0,IFLAG,IBIT )
5317312052 Jean*3131
1662f365b2 Andr*3132
3133
3134
13fa5cb63c Jean*3135
1662f365b2 Andr*3136 IF(IFLAG.GE.3) GO TO 75
13fa5cb63c Jean*3137
1662f365b2 Andr*3138 DO 9006 I = 1,IBIT
3139 ARG1(I) = 1. - X1(I)
3140 IF ( Z(I) .LT. 0.013 ) ARG1(I) =
3141 1 Z(I) * ( 0.25 - 0.09375 * Z(I) )
13fa5cb63c Jean*3142
1662f365b2 Andr*3143 ARG0(I) = 1. - X0(I)
3144 IF ( Z0(I) .LT. 0.013 ) ARG0(I) =
3145 1 Z0(I) * ( 0.25 - 0.09375 * Z0(I) )
13fa5cb63c Jean*3146
1662f365b2 Andr*3147 ARG1(I) = ARG1(I) * ( 1.+X0(I) )
3148 ARG0(I) = ARG0(I) * ( 1.+X1(I) )
3149 DX(I) = X1(I) - X0(I)
3150 ARG1(I) = ARG1(I) / ARG0(I)
3151 ARG0(I) = -DX(I) / ( 1. + X1(I)*X0(I) )
3152 ARG0(I) = ATAN( ARG0(I) )
3153 ARG1(I) = LOG( ARG1(I) )
3154 PSI2(I) = 2. * ARG0(I) + ARG1(I)
3155 PSI2(I) = PSI2(I) + DX(I)
3156 9006 CONTINUE
13fa5cb63c Jean*3157
471f7f004b Andr*3158 indx = 0
1662f365b2 Andr*3159 DO 9008 I = 1,IRUN
3160 IF( INTSTB(I).EQ.1 ) THEN
471f7f004b Andr*3161 indx = indx + 1
3162 VPSIM(I) = PSI2(indx)
3163 VX(I) = X1(indx)
3164 VXS(I) = X0(indx)
1662f365b2 Andr*3165 ENDIF
3166 9008 CONTINUE
13fa5cb63c Jean*3167
1662f365b2 Andr*3168
3169
3170
13fa5cb63c Jean*3171
1662f365b2 Andr*3172 IF(IFLAG.EQ.2) GO TO 100
13fa5cb63c Jean*3173
1662f365b2 Andr*3174 75 CONTINUE
3175 DO 9010 I = 1,IBIT
3176 ARG1(I) = 1. - Y1(I)
3177 IF( Z(I) .LT. 0.0065 ) ARG1(I) =
3178 1 Z(I) * ( 0.5 - 0.625 * Z(I) )
13fa5cb63c Jean*3179
1662f365b2 Andr*3180 ARG0(I) = 1. - Y0(I)
3181 IF( Z0(I) .LT. 0.0065 ) ARG0(I) =
3182 1 Z0(I) * ( 0.5 - 0.625 * Z0(I) )
13fa5cb63c Jean*3183
1662f365b2 Andr*3184 ARG1(I) = ARG1(I) * ( 1. + Y0(I) )
3185 ARG0(I) = ARG0(I) * ( 1. + Y1(I) )
3186 ARG1(I) = ARG1(I) / ARG0(I)
3187 PSI2(I) = LOG( ARG1(I) )
3188 PSI2(I) = PSI2(I) - Y1(I) + Y0(I)
3189 9010 CONTINUE
13fa5cb63c Jean*3190
471f7f004b Andr*3191 indx = 0
1662f365b2 Andr*3192 DO 9012 I = 1,IRUN
3193 IF( INTSTB(I).EQ.1 ) THEN
471f7f004b Andr*3194 indx = indx + 1
3195 VPSIH(I) = PSI2(indx)
3196 VY(I) = Y1(indx)
3197 VYS(I) = Y0(indx)
1662f365b2 Andr*3198 ENDIF
3199 9012 CONTINUE
13fa5cb63c Jean*3200
1662f365b2 Andr*3201
3202
3203
13fa5cb63c Jean*3204
1662f365b2 Andr*3205 100 CONTINUE
3206 IBIT = 0
3207 DO 9114 I = 1,IRUN
3208 IF(VZZ(I).GT.-1.e-7)THEN
3209 IBIT = IBIT + 1
3210 ENDIF
3211 9114 CONTINUE
3212 DO 9014 I = 1,IRUN
3213 IF(VZZ(I).GT.-1.e-7)THEN
3214 INTSTB(I) = 1
3215 ELSE
3216 INTSTB(I) = 0
3217 ENDIF
3218 9014 CONTINUE
3219 IF(IBIT.LE.0) GO TO 300
471f7f004b Andr*3220 indx = 0
06d4643e1f Jean*3221 #ifdef FIZHI_CRAY
1662f365b2 Andr*3222
3223 #endif
3224 DO 9016 I = 1,IRUN
3225 IF (INTSTB(I).EQ.1)THEN
471f7f004b Andr*3226 indx = indx + 1
3227 Z(indx) = VZZ(I)
3228 Z0(indx) = ZZ0(I)
3229 ARG1(indx) = VZH(I)
1662f365b2 Andr*3230 ENDIF
3231 9016 CONTINUE
06d4643e1f Jean*3232 #ifdef FIZHI_CRAY
1662f365b2 Andr*3233
3234 #endif
3235
3236 DO 9018 I = 1,IBIT
3237 HZ(I) = 1. / ARG1(I)
3238 Z1(I) = Z(I)
3239 Z2(I) = ZWM
13fa5cb63c Jean*3240
1662f365b2 Andr*3241 IF ( Z(I) .GT. ZWM ) THEN
3242 Z1(I) = ZWM
3243 Z2(I) = Z(I)
3244 ENDIF
13fa5cb63c Jean*3245
1662f365b2 Andr*3246 IF ( Z0(I) .GT. Z0M ) THEN
3247 Z0(I) = Z0M
3248 INTZ0(I) = 1
3249 ELSE
3250 INTZ0(I) = 0
3251 ENDIF
13fa5cb63c Jean*3252
1662f365b2 Andr*3253 X1NUM(I) = 1. + 5. * Z1(I)
3254 X0NUM(I) = 1. + 5. * Z0(I)
3255 X1DEN(I) = 1. / (1. + CM1 * (X1NUM(I) * Z1(I)) )
3256 X0DEN(I) = 1. + CM1 * (X0NUM(I) * Z0(I))
13fa5cb63c Jean*3257
1662f365b2 Andr*3258 IF ( (INTZ0(I).EQ.1) .OR. (Z(I).GT.ZWM) )
3259 1 HZ(I) = Z1(I) / Z0(I)
3260 ARG1(I) = HZ(I)*HZ(I)*X0DEN(I)*X1DEN(I)
3261 ARG1(I) = LOG( ARG1(I) )
3262 ARG1(I) = 0.5 * ARG1(I)
3263 ARG0(I) = (Z1(I) + 0.1) * (Z0(I) + 0.1)
3264 ARG0(I) = CM3 + ARG0(I) * CM4
3265 ARG0(I) = ( Z1(I) - Z0(I) ) / ARG0(I)
3266 ARG0(I) = ATAN( ARG0(I) )
3267 TEMP(I) = ARG1(I) + CM5 * ARG0(I)
13fa5cb63c Jean*3268
1662f365b2 Andr*3269 X0(I) = X0NUM(I) / X0DEN(I)
3270 IF ( INTZ0(I).EQ.1 ) X0(I) = 0.
3271 Z2ZWM(I) = Z2(I) * RZWM
3272 9018 CONTINUE
13fa5cb63c Jean*3273
1662f365b2 Andr*3274
3275
3276
13fa5cb63c Jean*3277
1662f365b2 Andr*3278 IF( IFLAG.GE.3 ) GO TO 225
13fa5cb63c Jean*3279
1662f365b2 Andr*3280 DO 9020 I = 1,IBIT
3281 X1(I) = X1NUM(I) * X1DEN(I)
3282 ARG1(I) = LOG( Z2ZWM(I) )
3283 PSI2(I) = TEMP(I) + CM6 * ARG1(I)
3284 9020 CONTINUE
13fa5cb63c Jean*3285
471f7f004b Andr*3286 indx = 0
1662f365b2 Andr*3287 DO 9030 I = 1,IRUN
3288 IF( INTSTB(I).EQ.1 ) THEN
471f7f004b Andr*3289 indx = indx + 1
3290 VPSIM(I) = PSI2(indx)
3291 VX(I) = X1(indx)
3292 VXS(I) = X0(indx)
1662f365b2 Andr*3293 ENDIF
3294 9030 CONTINUE
13fa5cb63c Jean*3295
1662f365b2 Andr*3296
3297
3298
13fa5cb63c Jean*3299
1662f365b2 Andr*3300 IF(IFLAG.EQ.2)GO TO 300
13fa5cb63c Jean*3301
1662f365b2 Andr*3302 225 CONTINUE
3303 DO 9024 I = 1,IBIT
3304 Y1DEN(I) = 1. + CM1 * ( X1NUM(I) * Z(I) )
3305 Y1(I) = X1NUM(I) / Y1DEN(I)
3306 ARG1(I) = CM7 * Z2ZWM(I) / ( CM2 + Z2(I) )
3307 ARG0(I) = 6.
3308 IF ( Z2(I) .GT. ZCM ) THEN
3309 Y1(I) = YCM
3310 ARG1(I) = Z2(I) * RZCM
3311 ARG0(I) = YCM
3312 TEMP(I) = TEMP(I) + CM8
3313 ENDIF
3314 ARG1(I) = LOG( ARG1(I) )
3315 PSI2(I) = TEMP(I) + ARG0(I) * ARG1(I)
3316 9024 CONTINUE
13fa5cb63c Jean*3317
471f7f004b Andr*3318 indx = 0
1662f365b2 Andr*3319 DO 9026 I = 1,IRUN
3320 IF( INTSTB(I).EQ.1 ) THEN
471f7f004b Andr*3321 indx = indx + 1
3322 VPSIH(I) = PSI2(indx)
3323 VY(I) = Y1(indx)
3324 VYS(I) = X0(indx)
1662f365b2 Andr*3325 ENDIF
3326 9026 CONTINUE
13fa5cb63c Jean*3327
1662f365b2 Andr*3328 300 CONTINUE
13fa5cb63c Jean*3329
1662f365b2 Andr*3330 RETURN
3331 END
13fa5cb63c Jean*3332
3333
3334
1662f365b2 Andr*3335 SUBROUTINE TRBLEN (STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,
3336 1 QXLM,NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
3337 2 DZITRP,STBFCN,XL0,Q1,WRKIT1,WRKIT2,WRKIT3,WRKIT4,irun)
3338
3339
3340
3341
3342
3343
3344
3345
3346
3347
3348
3349
3350
3351
3352
3353
3354
3355
3356
3357
3358
3359
3360
3361
3362
3363
3364
3365
3366
3367
453e05dab3 Andr*3368 implicit none
3369
3370
3371 integer irun,nlev,init,lmin,lminq,lminq1
a456aa407c Andr*3372 _RL cp
3373 _RL STRT(irun,NLEV),DW2(irun,NLEV),DZ3(irun,NLEV)
3374 _RL Q(irun,NLEV),VKZM(irun,NLEV-1),VKZE(irun,NLEV-1)
3375 _RL DTHV(irun,NLEV),DPK(irun,NLEV),DU(irun,NLEV)
3376 _RL DV(irun,NLEV)
3377 _RL QXLM(irun,NLEV-1),XL(irun,NLEV-1)
3378 _RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
3379 _RL XL0(irun,nlev),Q1(irun,nlev-1)
3380 _RL WRKIT1(irun,nlev-1),WRKIT2(irun,nlev-1)
3381 _RL WRKIT3(irun,nlev-1)
3382 _RL WRKIT4(irun,nlev-1)
453e05dab3 Andr*3383 INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
3384
3385
a456aa407c Andr*3386 _RL rf1,rf2,e5,d4,d1,rfc,ricr,alpha,dzcnv,xl0cnv,xl0min
3387 _RL clmt,clmt53
1662f365b2 Andr*3388 PARAMETER ( RF1 = 0.2340678 )
3389 PARAMETER ( RF2 = 0.2231172 )
5317312052 Jean*3390 PARAMETER ( E5 = 49.66 )
1662f365b2 Andr*3391 PARAMETER ( D4 = 2.6532122E-2 )
3392 PARAMETER ( D1 = D4 * E5 )
3393 PARAMETER ( RFC = 0.1912323 )
3394 PARAMETER ( RICR = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 ) )
3395 PARAMETER ( ALPHA = 0.1 )
3396 PARAMETER ( DZCNV = 100. )
3397 PARAMETER ( XL0CNV = DZCNV * ALPHA )
3398 PARAMETER ( XL0MIN = 1. )
5317312052 Jean*3399 PARAMETER ( CLMT = 0.23 )
1662f365b2 Andr*3400 PARAMETER ( CLMT53 = 5. * CLMT / 3. )
453e05dab3 Andr*3401
3402 integer ibit,nlevm1,nlevp1,istnlv,istnm1,nlevml,istnml,Lp1
3403 integer istnmq,istlmq,lminp,lm1,lmin1
3404 integer i,L,LL
13fa5cb63c Jean*3405
1662f365b2 Andr*3406 NLEVM1 = NLEV - 1
3407 NLEVP1 = NLEV + 1
3408 ISTNLV = irun * NLEV
3409 ISTNM1 = irun * NLEVM1
13fa5cb63c Jean*3410
1662f365b2 Andr*3411 IF ( INIT.EQ.2 ) GO TO 1200
13fa5cb63c Jean*3412
1662f365b2 Andr*3413
3414
5317312052 Jean*3415 DO L =1,NLEV
3416 DO I =1,irun
3417 STBFCN(I,L) = STRT(I,L) - RICR * DW2(I,L)
3418 INT1(I,L) = 0
3419 IF( STBFCN(I,L).LE.0. ) INT1(I,L) = 1
3420 ENDDO
3421 ENDDO
3422
3423 DO L =1,NLEVM1
3424 DO I =1,irun
3425 INT2(I,L) = 0
3426 IF( (INT1(I,L).EQ.1) .NEQV. (INT1(I,L+1).EQ.1) ) INT2(I,L) = 1
3427 ENDDO
3428 ENDDO
13fa5cb63c Jean*3429
1662f365b2 Andr*3430 DO 40 LMIN = 1,NLEV
3431 IBIT = 0
3432 DO 30 I=1,irun
3433 IBIT = IBIT + INT1(I,LMIN)
3434 30 CONTINUE
3435 IF(IBIT.GE.1) GO TO 50
3436 40 CONTINUE
3437 LMIN = NLEVP1
3438 50 CONTINUE
3439 LMIN = 1
13fa5cb63c Jean*3440
5317312052 Jean*3441 DO L =1,NLEVM1
3442 DO I =1,irun
3443 XL0(I,L) = 0.
3444 ENDDO
3445 ENDDO
1662f365b2 Andr*3446 DO 70 I=1,irun
3447 XL0(I,NLEV) = DZ3(I,NLEV)
3448 70 CONTINUE
13fa5cb63c Jean*3449
1662f365b2 Andr*3450 IF(LMIN.GE.NLEVP1) GOTO 1100
3451 LMIN1 = LMIN - 1
3452 IF(LMIN1.EQ.0) LMIN1 = 1
3453 NLEVML = NLEV - LMIN1
3454 ISTNML = irun*NLEVML
3455 CALL TRBITP ( STBFCN(1,LMIN1),INT2(1,LMIN1),DTHV(1,LMIN1),
13fa5cb63c Jean*3456 & DPK(1,LMIN1), DU(1,LMIN1), DV(1,LMIN1),
3457 & DZITRP(1,LMIN1), NLEVML,
3458 & WRKIT1,WRKIT2,WRKIT3,WRKIT4,CP,irun )
1662f365b2 Andr*3459 LP1 = LMIN1 + 1
13fa5cb63c Jean*3460
5317312052 Jean*3461 DO L =LMIN1,NLEVM1
3462 DO I =1,irun
3463 INT2(I,L) = 0
3464 IF ( INT1(I,L).EQ.1 .OR. INT1(I,L+1).EQ.1 ) INT2(I,L) = 1
3465 IF ( INT2(I,L).EQ.1 )
13fa5cb63c Jean*3466 & XL0(I,L) = (0.5+DZITRP(I,L)) * DZ3(I,L+1)
5317312052 Jean*3467 ENDDO
3468 ENDDO
1662f365b2 Andr*3469 DO 90 I=1,irun
3470 INT2(I,NLEVM1) = INT1(I,NLEV)
3471 90 CONTINUE
13fa5cb63c Jean*3472
5317312052 Jean*3473 DO L =LMIN1,NLEVM1
3474 DO I =1,irun
3475 IF ( INT2(I,L).EQ.1 ) THEN
3476 XL0(I,L+1) = XL0(I,L+1) + ( (0.5-DZITRP(I,L)) * DZ3(I,L+1) )
3477 ENDIF
3478 ENDDO
3479 ENDDO
1662f365b2 Andr*3480 IF (LMIN.GT.1) GOTO 400
3481 DO 110 I=1,irun
3482 IF( INT1(I,1).EQ.1 ) XL0(I,1) = XL0(I,1) + DZ3(I,1)
3483 110 CONTINUE
3484 400 CONTINUE
13fa5cb63c Jean*3485
1662f365b2 Andr*3486 LMINP = LMIN + 1
3487 IF(LMINP.GT.NLEVM1) GOTO 550
3488 DO 500 L = LMINP,NLEVM1
3489 LM1 = L-1
3490 DO 120 I = 1,irun
3491 IF( INT1(I,LM1).EQ.1 ) XL0(I,L) = XL0(I,L) + XL0(I,LM1)
3492 120 CONTINUE
3493 500 CONTINUE
3494 550 CONTINUE
3495 IF(LMIN.GT.NLEVM1) GOTO 600
3496 DO 130 I = 1,irun
3497 IF( INT1(I,NLEVM1).EQ.1 .AND. INT1(I,NLEV).EQ.1 ) THEN
3498 XL0(I,NLEV) = XL0(I,NLEV) + XL0(I,NLEVM1)
3499 ENDIF
3500 130 CONTINUE
3501 IF(LMIN.GT.NLEV) GOTO 1100
3502 600 CONTINUE
3503 DO 1000 LL = LMIN,NLEV-1
3504 L = NLEVM1 + LMIN - LL
3505 LP1 = L+1
3506 DO 140 I = 1,irun
3507 IF( INT1(I,LP1).EQ.1 ) THEN
3508 IF( INT1(I,L) .EQ.1 ) THEN
3509 XL0(I,L) = XL0(I,LP1)
3510 ELSE
3511 XL0(I,L) = XL0(I,L) + XL0(I,LP1)
3512 ENDIF
3513 ENDIF
3514 140 CONTINUE
3515 1000 CONTINUE
3516 1100 CONTINUE
13fa5cb63c Jean*3517
5317312052 Jean*3518 DO L =1,NLEV
3519 DO I =1,irun
3520 IF ( XL0(I,L).LT.XL0CNV ) XL0(I,L) = XL0CNV
3521 ENDDO
3522 ENDDO
13fa5cb63c Jean*3523
1662f365b2 Andr*3524
3525
3526
13fa5cb63c Jean*3527
1662f365b2 Andr*3528 IF(INIT.EQ.1) GOTO 1400
13fa5cb63c Jean*3529
1662f365b2 Andr*3530 IF(LMINQ.GT.1) THEN
3531 ISTLMQ = irun * LMINQ1
5317312052 Jean*3532 DO L =1,LMINQ1
3533 DO I =1,irun
3534 INT2(I,L) = 1 - INT1(I,L)
3535 ENDDO
3536 ENDDO
1662f365b2 Andr*3537 ENDIF
3538 IF(LMINQ.LT.NLEV) THEN
3539 ISTNMQ = irun * (NLEV-LMINQ)
5317312052 Jean*3540 DO L =LMINQ,NLEVM1
3541 DO I =1,irun
3542 IF( INT1(I,L).EQ.0 ) THEN
3543 XL0(I,L) = Q(I,L) / XL0(I,L)
3544 XL0(I,L) = XL0(I,L) * XL0(I,L) + 1.0E-20
3545 XL0(I,L) = STBFCN(I,L) + XL0(I,L)
3546 XL0(I,L) = SQRT( XL0(I,L) )
3547 XL0(I,L) = Q(I,L) / XL0(I,L)
3548 ENDIF
3549 INT2(I,L) = 0
3550 IF( XL0(I,L).LT.XL0MIN ) INT2(I,L) = 1
3551 ENDDO
3552 ENDDO
1662f365b2 Andr*3553 ENDIF
13fa5cb63c Jean*3554
1662f365b2 Andr*3555 1200 CONTINUE
13fa5cb63c Jean*3556
1662f365b2 Andr*3557 IF(INIT.EQ.2) THEN
5317312052 Jean*3558 DO L =1,NLEVM1
3559 DO I =1,irun
3560 INT2(I,L) = 1 - INT1(I,L)
3561 IF ( INT2(I,L).EQ.1 ) XL0(I,L) = XL0MIN
3562 ENDDO
3563 ENDDO
1662f365b2 Andr*3564 ENDIF
5317312052 Jean*3565 DO L =1,NLEVM1
3566 DO I =1,irun
3567 IF ( INT2(I,L).EQ.1 ) XL0(I,L) = XL0MIN
3568 ENDDO
3569 ENDDO
13fa5cb63c Jean*3570
1662f365b2 Andr*3571
3572
3573
13fa5cb63c Jean*3574
1662f365b2 Andr*3575 1400 CONTINUE
13fa5cb63c Jean*3576
5317312052 Jean*3577 DO L =1,NLEVM1
3578 DO I =1,irun
3579 XL(I,L) = XL0(I,L) * VKZE(I,L) / ( XL0(I,L)+VKZE(I,L) )
3580 ENDDO
3581 ENDDO
13fa5cb63c Jean*3582
1662f365b2 Andr*3583
3584
3585
13fa5cb63c Jean*3586
1662f365b2 Andr*3587 IF(INIT.EQ.1) GOTO 1700
3588 ISTNMQ = irun * (NLEV-LMINQ1)
5317312052 Jean*3589 DO L =LMINQ1,NLEVM1
3590 DO I =1,irun
3591 Q1(I,L) = Q(I,L)
3592 INT1(I,L) = 0
3593 IF ( Q(I,L).LE.Q(I,L+1) ) INT1(I,L) = 1
3594 IF ( INT1(I,L).EQ.1 ) THEN
3595 XL0(I,L) = XL0(I,L+1)
3596 Q1(I,L) = Q(I,L+1)
3597 ENDIF
3598 ENDDO
3599 ENDDO
13fa5cb63c Jean*3600
5317312052 Jean*3601 DO L =LMINQ1,NLEVM1
3602 DO I =1,irun
3603 QXLM(I,L) = XL0(I,L)*VKZM(I,L)
13fa5cb63c Jean*3604 & / ( XL0(I,L)+VKZM(I,L) )
5317312052 Jean*3605 QXLM(I,L) = CLMT53 * Q1(I,L)*QXLM(I,L)
3606 ENDDO
3607 ENDDO
13fa5cb63c Jean*3608
1662f365b2 Andr*3609 1700 CONTINUE
13fa5cb63c Jean*3610
1662f365b2 Andr*3611 RETURN
3612 END
13fa5cb63c Jean*3613
3614
3615
1662f365b2 Andr*3616 SUBROUTINE TRBITP ( STBFCN,INTCHG,DTHV,DPK,DU,DV,DZITRP,NLEV,
13fa5cb63c Jean*3617 & AAA,BBB,CCC,DDD,CP,irun )
1662f365b2 Andr*3618
3619
3620
3621
3622
3623
3624
3625
3626
3627
3628
3629
3630
3631
3632
3633
3634
453e05dab3 Andr*3635 implicit none
3636
3637
3638 integer irun,nlev
a456aa407c Andr*3639 _RL cp
3640 _RL STBFCN(irun,NLEV+1)
453e05dab3 Andr*3641 integer INTCHG(irun,NLEV)
a456aa407c Andr*3642 _RL DTHV(irun,NLEV+1),DPK(irun,NLEV+1)
3643 _RL DU(irun,NLEV+1),DV(irun,NLEV+1)
3644 _RL DZITRP(irun,NLEV+1)
3645 _RL AAA(irun,NLEV),BBB(irun,NLEV)
3646 _RL CCC(irun,NLEV),DDD(irun,NLEV)
453e05dab3 Andr*3647
3648
a456aa407c Andr*3649 _RL rf1,rf2,e5,d4,d1,rfc,ricr
1662f365b2 Andr*3650 PARAMETER ( RF1 = 0.2340678 )
3651 PARAMETER ( RF2 = 0.2231172 )
5317312052 Jean*3652 PARAMETER ( E5 = 49.66 )
1662f365b2 Andr*3653 PARAMETER ( D4 = 2.6532122E-2 )
3654 PARAMETER ( D1 = D4 * E5 )
3655 PARAMETER ( RFC = 0.1912323 )
3656 PARAMETER ( RICR = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 ) )
3657
453e05dab3 Andr*3658 integer istnlv
5317312052 Jean*3659 integer I,L
13fa5cb63c Jean*3660
1662f365b2 Andr*3661
3662
3663
3664
13fa5cb63c Jean*3665
1662f365b2 Andr*3666 ISTNLV = irun*NLEV
5317312052 Jean*3667 DO L =1,NLEV
3668 DO I =1,irun
3669 DZITRP(I,L) = 0.
3670 ENDDO
3671 ENDDO
3672 DO L =1,NLEV
3673 DO I =1,irun
3674 IF ( INTCHG(I,L).EQ.1 ) THEN
3675 DDD(I,L) = ( CP *(DTHV(I,L+1)*DPK(I,L)
13fa5cb63c Jean*3676 & + DTHV(I,L)*DPK(I,L+1)) )
3677 & - ( (2.*RICR) * ( DU(I,L+1)* DU(I,L)
3678 & + DV(I,L+1)* DV(I,L)) )
5317312052 Jean*3679 AAA(I,L) = STBFCN(I,L) + STBFCN(I,L+1)
3680 BBB(I,L) = STBFCN(I,L) - STBFCN(I,L+1)
3681 CCC(I,L) = 1. / BBB(I,L)
3682 DZITRP(I,L) = AAA(I,L) * CCC(I,L)
3683 AAA(I,L) = AAA(I,L) - DDD(I,L)
3684 DDD(I,L) = ( DDD(I,L) * DDD(I,L) )
13fa5cb63c Jean*3685 & - 4. * (STBFCN(I,L+1) * STBFCN(I,L) )
5317312052 Jean*3686 DDD(I,L) = DDD(I,L)*CCC(I,L)*CCC(I,L)
3687 DDD(I,L) = SQRT( DDD(I,L) )
3688 ENDIF
13fa5cb63c Jean*3689
5317312052 Jean*3690 IF( INTCHG(I,L).EQ.1 .AND. AAA(I,L).NE.0. ) THEN
3691 DZITRP(I,L) = ( BBB(I,L)*(1.-DDD(I,L)) ) / AAA(I,L)
3692 ENDIF
13fa5cb63c Jean*3693
5317312052 Jean*3694 DZITRP(I,L) = 0.5 * DZITRP(I,L)
3695 ENDDO
3696 ENDDO
13fa5cb63c Jean*3697
1662f365b2 Andr*3698 RETURN
3699 END
13fa5cb63c Jean*3700
3701
3702
1662f365b2 Andr*3703 SUBROUTINE TRBL20 (RI,STRT,DW2,XL,ZKM,ZKH,QE,QQE,INTSTB,NLEV,
3704 1 nlay,irun)
3705
3706
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717
3718
3719
3720
3721
3722
3723
3724
453e05dab3 Andr*3725 implicit none
3726
3727
3728 integer nlev,nlay,irun
5317312052 Jean*3729 _RL RI(irun*NLEV,1),STRT(irun*NLEV,1),DW2(irun*NLEV,1)
3730 _RL XL(irun*NLEV,1),ZKM(irun*NLEV,1), ZKH(irun*NLEV,1)
3731 _RL QE(irun*NLEV,1),QQE(irun*NLEV,1)
3732 INTEGER INTSTB(irun*NLEV,1)
3733 _RL EE(irun*(nlay-1),1),RF(irun*(nlay-1),1)
453e05dab3 Andr*3734
3735
a456aa407c Andr*3736 _RL b1,b2,d3,rf1,rf2,d3b2,d2,e5,d4,d1,d1half,d2half
5317312052 Jean*3737 _RL rfc,ricr,ch,cm
3738 PARAMETER ( B1 = 16.6 )
1662f365b2 Andr*3739 PARAMETER ( B2 = 10.1 )
3740 PARAMETER ( D3 = 0.29397643 )
3741 PARAMETER ( RF1 = 0.2340678 )
3742 PARAMETER ( RF2 = 0.2231172 )
3743 PARAMETER ( D3B2 = D3 / RF1 )
3744 PARAMETER ( D2 = RF1 )
5317312052 Jean*3745 PARAMETER ( E5 = 49.66 )
1662f365b2 Andr*3746 PARAMETER ( D4 = 2.6532122E-2 )
3747 PARAMETER ( D1 = D4 * E5 )
3748 PARAMETER ( D1HALF = 0.5 * D1 )
3749 PARAMETER ( D2HALF = 0.5 * D2 )
3750 PARAMETER ( RFC = 0.1912323 )
3751 PARAMETER ( RICR = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 ) )
3752 PARAMETER ( CH = 2.5828674 )
3753 PARAMETER ( CM = CH / D1 )
3754
453e05dab3 Andr*3755 integer istnlv
3756 integer i
3757
1662f365b2 Andr*3758 ISTNLV = irun * NLEV
13fa5cb63c Jean*3759
1662f365b2 Andr*3760
3761
3762
13fa5cb63c Jean*3763
1662f365b2 Andr*3764 DO 10 I=1,ISTNLV
3765 EE(I,1) = D1HALF * RI(I,1) + D2HALF
3766 RF(I,1) = EE(I,1)* EE(I,1)
3767 RF(I,1) = RF(I,1)- D3*RI(I,1)
3768 RF(I,1) = SQRT( RF(I,1) )
3769 RF(I,1) = EE(I,1) - RF(I,1)
13fa5cb63c Jean*3770
1662f365b2 Andr*3771 IF( RI(I,1).LE.1.e-4 .AND. RI(I,1).GE.-1.e-4 ) THEN
3772 RF(I,1) = D3B2*RI(I,1)
3773 ENDIF
13fa5cb63c Jean*3774
1662f365b2 Andr*3775
3776
3777
3778
13fa5cb63c Jean*3779
1662f365b2 Andr*3780 IF( RI(I,1).LT.RICR .AND. RF(I,1).LT.RFC ) THEN
3781 ZKH(I,1) = ( RFC-RF(I,1) ) / (1.-RF(I,1))
3782 ZKM(I,1) = CM * (RF1-RF(I,1))
3783 ZKM(I,1) = ZKH(I,1)*ZKM(I,1) / (RF2-RF(I,1))
3784 ZKH(I,1) = CH *ZKH(I,1)
3785 QE(I,1) = ZKM(I,1)*DW2(I,1) - ZKH(I,1)*STRT(I,1)
3786 ENDIF
13fa5cb63c Jean*3787
1662f365b2 Andr*3788 IF( QE(I,1).LT.1.e-14 ) THEN
3789 INTSTB(I,1) = 0
3790 QE(I,1) = 0.
3791 ELSE
3792 INTSTB(I,1) = 1
3793 QE(I,1) = B1*QE(I,1)
3794 QE(I,1) = SQRT( QE(I,1) )
3795 QE(I,1) = XL(I,1)*QE(I,1)
3796 ENDIF
3797 QQE(I,1) = 0.5 * QE(I,1) * QE(I,1)
3798 10 CONTINUE
13fa5cb63c Jean*3799
1662f365b2 Andr*3800 RETURN
3801 END
13fa5cb63c Jean*3802
3803
3804
1662f365b2 Andr*3805 SUBROUTINE TRBL25(Q,XL,STRT,DW2,INTSTB,INTQ,ZKM,ZKH,P3,NLEV,
3806 1 nlay,irun)
3807
3808
3809
3810
3811
3812
3813
3814
3815
3816
3817
3818
3819
3820
3821
3822
3823
3824
3825
453e05dab3 Andr*3826 implicit none
3827
3828
3829 integer nlev,nlay,irun
5317312052 Jean*3830 _RL Q(irun*NLEV,1),XL(irun*NLEV,1),STRT(irun*NLEV,1)
3831 _RL DW2(irun*NLEV,1)
3832 INTEGER INTSTB(irun*nlay,1), INTQ(irun*nlay,1)
3833 _RL ZKM(irun*NLEV,1),ZKH(irun*NLEV,1),P3(irun*NLEV,1)
453e05dab3 Andr*3834
3835
a456aa407c Andr*3836 _RL a1,a2,a4,c1,a5,a3,b1,b2,b3,ff2,ff3,ff4
1662f365b2 Andr*3837 PARAMETER ( A1 = 0.92 )
3838 PARAMETER ( A2 = 0.74 )
3839 PARAMETER ( A4 = 6. * A1 * A1)
3840 PARAMETER ( C1 = 0.08 )
3841 PARAMETER ( A5 = 3.*C1*(-1.) )
3842 PARAMETER ( A3 = A4 * A5*(-1.) )
5317312052 Jean*3843 PARAMETER ( B1 = 16.6 )
1662f365b2 Andr*3844 PARAMETER ( B2 = 10.1 )
5317312052 Jean*3845 PARAMETER ( B3 = 1. / B1 )
1662f365b2 Andr*3846 PARAMETER ( FF2 = 9. * A1 * A2 )
3847 PARAMETER ( FF3 = (3.*A2*B2) - (9.*A2*A2 ) )
3848 PARAMETER ( FF4 = (3.*A2*B2) + (12.*A1*A2 ) )
3849
5317312052 Jean*3850 _RL F2(irun*(nlay-1),1),F3(irun*(nlay-1),1)
3851 _RL F4(irun*(nlay-1),1),XQ(irun*(nlay-1),1)
453e05dab3 Andr*3852
3853 integer istnlv
3854 integer i
13fa5cb63c Jean*3855
1662f365b2 Andr*3856 ISTNLV = irun * NLEV
13fa5cb63c Jean*3857
1662f365b2 Andr*3858
3859
3860
3861
13fa5cb63c Jean*3862
1662f365b2 Andr*3863 DO 10 I=1,ISTNLV
3864 IF( INTQ(I,1).EQ.1 .AND. INTSTB(I,1).EQ.0 ) THEN
3865 XQ(I,1) = XL(I,1) / Q(I,1)
3866 XQ(I,1) = XQ(I,1) * XQ(I,1)
3867 STRT(I,1) = XQ(I,1) * STRT(I,1)
3868 DW2(I,1) = XQ(I,1) * DW2(I,1)
3869 F2(I,1) = 1.+FF2 * STRT(I,1)
3870 F3(I,1) = 1.+FF3 * STRT(I,1)
3871 F4(I,1) = 1.+FF4 * STRT(I,1)
3872 ZKH(I,1) = (F4(I,1) * F2(I,1))
13fa5cb63c Jean*3873 & + A4 * (F3(I,1) * DW2(I,1))
1662f365b2 Andr*3874 ZKH(I,1) = (F2(I,1) + A3*DW2(I,1))
13fa5cb63c Jean*3875 & / ZKH(I,1)
1662f365b2 Andr*3876 ZKM(I,1) = A1 * (F3(I,1)*ZKH(I,1)+A5)
13fa5cb63c Jean*3877 & / F2(I,1)
1662f365b2 Andr*3878 ZKH(I,1) = A2 * ZKH(I,1)
3879 P3(I,1) = ZKH(I,1)*STRT(I,1) + B3
3880 P3(I,1) = 2. * ( ZKM(I,1)*DW2(I,1) - P3(I,1) )
3881 P3(I,1) = P3(I,1)*Q(I,1)
13fa5cb63c Jean*3882
1662f365b2 Andr*3883 ENDIF
3884 10 CONTINUE
13fa5cb63c Jean*3885
1662f365b2 Andr*3886 RETURN
3887 END
13fa5cb63c Jean*3888
3889
3890
1662f365b2 Andr*3891 SUBROUTINE TRBDIF ( XX1,XX2,RHOKDZ,FLXFAC,DXX1G,DXX2G,NLEV,
13fa5cb63c Jean*3892 & ITYPE,EPSL,irun )
3893
1662f365b2 Andr*3894
3895
3896
3897
3898
3899
3900
3901
3902
3903
3904
3905
3906
3907
3908
3909
3910
3911
3912
3913
3914
3915
3916
3917
3918
3919
3920
3921
3922
3923
453e05dab3 Andr*3924 implicit none
3925
3926
3927 integer nlev,itype,irun
a456aa407c Andr*3928 _RL XX1(irun,NLEV+1),XX2(irun,NLEV+1)
3929 _RL RHOKDZ(irun,NLEV),FLXFAC(irun,NLEV+1)
3930 _RL DXX1G(irun),DXX2G(irun)
3931 _RL epsl
13fa5cb63c Jean*3932
a456aa407c Andr*3933 _RL AA(irun,nlev), BB(irun,nlev), CC(irun,nlev+1)
5317312052 Jean*3934 integer NLEVP1,NLEVX
3935 integer I,L
13fa5cb63c Jean*3936
1662f365b2 Andr*3937 NLEVP1 = NLEV + 1
5317312052 Jean*3938 NLEVX = NLEV - 1
3939 IF(ITYPE.EQ.2) NLEVX = NLEV
13fa5cb63c Jean*3940
1662f365b2 Andr*3941
13fa5cb63c Jean*3942
1662f365b2 Andr*3943 DO 10 I=1,irun
3944 CC(I,1) = 0.
3945 10 CONTINUE
5317312052 Jean*3946 DO L =1,NLEVX
3947 DO I =1,irun
3948 CC(I,L+1) = RHOKDZ(I,L) * FLXFAC(I,L+1)
3949 ENDDO
3950 ENDDO
3951 DO L =1,NLEV
3952 DO I =1,irun
3953 BB(I,L) = RHOKDZ(I,L) * FLXFAC(I,L)
3954 AA(I,L) = 1. + CC(I,L) + BB(I,L)
3955 ENDDO
3956 ENDDO
13fa5cb63c Jean*3957
1662f365b2 Andr*3958
3959 IF(ITYPE.EQ.1) THEN
5317312052 Jean*3960 DO L =1,NLEV
3961 DO I =1,irun
3962 AA(I,L) = AA(I,L) - XX2(I,L)
3963 ENDDO
3964 ENDDO
1662f365b2 Andr*3965 ENDIF
13fa5cb63c Jean*3966
1662f365b2 Andr*3967
3968 CALL VTRI0(AA,BB,CC,XX1,XX1,NLEV,irun)
13fa5cb63c Jean*3969
1662f365b2 Andr*3970 IF(ITYPE.EQ.2) THEN
3971
13fa5cb63c Jean*3972
1662f365b2 Andr*3973 DO 50 I=1,irun
3974 DXX1G(I) = CC(I,NLEVP1) * ( XX1(I,NLEV)-XX1(I,NLEVP1) )
3975 50 CONTINUE
13fa5cb63c Jean*3976
1662f365b2 Andr*3977
3978 CALL VTRI1(AA,BB,XX2,NLEV,irun)
3979 DO 60 I=1,irun
3980 DXX2G(I) = CC(I,NLEVP1) * ( XX2(I,NLEV)-XX2(I,NLEVP1) )
3981 60 CONTINUE
3982 ENDIF
13fa5cb63c Jean*3983
1662f365b2 Andr*3984
13fa5cb63c Jean*3985
1662f365b2 Andr*3986 IF(ITYPE.EQ.3) CALL VTRI2 (AA,BB,CC,XX2,XX2,NLEV,irun)
13fa5cb63c Jean*3987
1662f365b2 Andr*3988
3989 IF(ITYPE.EQ.1) THEN
5317312052 Jean*3990 DO L =1,NLEV
3991 DO I =1,irun
3992 IF ( XX1(I,L).LT.EPSL ) XX1(I,L) = 0.
3993 ENDDO
3994 ENDDO
1662f365b2 Andr*3995 ENDIF
13fa5cb63c Jean*3996
1662f365b2 Andr*3997 RETURN
3998 END
13fa5cb63c Jean*3999
4000
4001
1662f365b2 Andr*4002 SUBROUTINE VTRI0 ( A,B,C,F,Y,K,irun)
453e05dab3 Andr*4003 implicit none
4004
4005 integer k,irun
a456aa407c Andr*4006 _RL A(irun,K),B(irun,K),C(irun,K),Y(irun,K+1)
4007 _RL F(irun,K)
453e05dab3 Andr*4008
4009 integer i,L,Lm1
13fa5cb63c Jean*4010
1662f365b2 Andr*4011 DO 9000 I = 1,irun
4012 A(I,1) = 1. / A(I,1)
4013 9000 CONTINUE
13fa5cb63c Jean*4014
1662f365b2 Andr*4015 DO 100 L = 2,K
4016 LM1 = L - 1
4017 DO 9002 I = 1,irun
4018 C(I,L) = C(I,L) * A(I,LM1)
4019 A(I,L) = 1. / ( A(I,L) - B(I,LM1) * C(I,L) )
4020 F(I,L) = F(I,L) + F(I,LM1) * C(I,L)
4021 9002 CONTINUE
4022 100 CONTINUE
13fa5cb63c Jean*4023
1662f365b2 Andr*4024 DO 200 L = K,1,-1
4025 DO 9004 I = 1,irun
4026 Y(I,L) = (F(I,L) + B(I,L) * Y(I,L+1)) * A(I,L)
4027 9004 CONTINUE
4028 200 CONTINUE
13fa5cb63c Jean*4029
1662f365b2 Andr*4030 RETURN
4031 END
13fa5cb63c Jean*4032
4033
4034
1662f365b2 Andr*4035 SUBROUTINE VTRI1 ( A,B,Y,K,irun)
453e05dab3 Andr*4036 implicit none
4037
4038 integer k,irun
a456aa407c Andr*4039 _RL A(irun,K),B(irun,K),Y(irun,K+1)
453e05dab3 Andr*4040
4041 integer i,L
13fa5cb63c Jean*4042
1662f365b2 Andr*4043 DO 200 L = K,1,-1
4044 DO 9000 I = 1,irun
4045 Y(I,L) = B(I,L) * Y(I,L+1) * A(I,L)
4046 9000 CONTINUE
4047 200 CONTINUE
13fa5cb63c Jean*4048
1662f365b2 Andr*4049 RETURN
4050 END
13fa5cb63c Jean*4051
4052
4053
1662f365b2 Andr*4054 SUBROUTINE VTRI2 ( A,B,C,F,Y,K,irun)
453e05dab3 Andr*4055 implicit none
4056
4057 integer k,irun
a456aa407c Andr*4058 _RL A(irun,K),B(irun,K),C(irun,K),F(irun,K)
4059 _RL Y(irun,K+1)
453e05dab3 Andr*4060
4061 integer i,L
13fa5cb63c Jean*4062
1662f365b2 Andr*4063 DO 100 L = 2,K
4064 DO 9000 I = 1,irun
4065 F(I,L) = F(I,L) + F(I,L-1) * C(I,L)
4066 9000 CONTINUE
4067 100 CONTINUE
13fa5cb63c Jean*4068
1662f365b2 Andr*4069 DO 200 L = K,1,-1
4070 DO 9002 I = 1,irun
4071 Y(I,L) = (F(I,L) + B(I,L) * Y(I,L+1)) * A(I,L)
4072 9002 CONTINUE
4073 200 CONTINUE
13fa5cb63c Jean*4074
1662f365b2 Andr*4075 RETURN
4076 END
13fa5cb63c Jean*4077
4078
4079
1662f365b2 Andr*4080 SUBROUTINE LINADJ ( NN,VRIB1,VRIB2,VWS1,VWS2,VZ1,VUSTAR,IWATER,
4081 1 VAPSIM, VAPSIHG,VPSIH,VPSIG,VX,VX0,VY,VY0,ITYPE,LWATER,IRUN,
4082 2 VDZETA,VDZ0,VDPSIM,VDPSIH,INTRIB,
4083 3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
4084 4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
4085 5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
13fa5cb63c Jean*4086
1662f365b2 Andr*4087
4088
4089
4090
4091
4092
4093
4094
4095
4096
4097
4098
4099
4100
4101
4102
4103
4104
4105
4106
4107
4108
4109
4110
4111
4112
4113
4114
4115
4116
4117
4118
4119
4120
453e05dab3 Andr*4121 implicit none
4122
4123
4124 integer nn,irun,itype
a456aa407c Andr*4125 _RL VRIB1(IRUN),VRIB2(IRUN)
4126 _RL VWS1(IRUN),VWS2(IRUN),VZ1(IRUN),VUSTAR(IRUN)
453e05dab3 Andr*4127 integer IWATER(IRUN)
a456aa407c Andr*4128 _RL VAPSIM(IRUN),VAPSIHG(IRUN)
4129 _RL VPSIH(IRUN),VPSIG(IRUN),VX(IRUN)
4130 _RL VX0(IRUN),VY(IRUN),VY0(IRUN)
453e05dab3 Andr*4131 LOGICAL LWATER
a456aa407c Andr*4132 _RL VDZETA(IRUN),VDZ0(IRUN),VDPSIM(IRUN)
4133 _RL VDPSIH(IRUN)
453e05dab3 Andr*4134 integer INTRIB(IRUN)
a456aa407c Andr*4135 _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
4136 _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
4137 _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
4138 _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
4139 _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
4140 _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
4141 _RL VDPSIMC(irun),VDPSIHC(irun)
453e05dab3 Andr*4142
4143
a456aa407c Andr*4144 _RL xx0max,prfac,xpfac,difsqt,ustz0s,h0byz0,usth0s
1662f365b2 Andr*4145 PARAMETER ( XX0MAX = 1.49821 )
4146 PARAMETER ( PRFAC = 0.595864 )
5317312052 Jean*4147 PARAMETER ( XPFAC = .55 )
1662f365b2 Andr*4148 PARAMETER ( DIFSQT = 3.872983E-3)
4149 PARAMETER ( USTZ0S = 0.2030325E-5)
4150 PARAMETER ( H0BYZ0 = 30.0 )
4151 PARAMETER ( USTH0S = H0BYZ0*USTZ0S )
4152
453e05dab3 Andr*4153 integer VINT1(irun),VINT2(irun)
a456aa407c Andr*4154 _RL getcon,vk,bmdl,b2uhs
453e05dab3 Andr*4155 integer i
13fa5cb63c Jean*4156
1662f365b2 Andr*4157 vk = getcon('VON KARMAN')
4158 BMDL = VK * XPFAC * PRFAC / DIFSQT
4159 B2UHS = BMDL * BMDL * USTH0S
4160
4161
4162
13fa5cb63c Jean*4163
1662f365b2 Andr*4164
4165 IF ( (ITYPE.EQ.1) .AND. LWATER ) THEN
4166 DO 9000 I = 1,IRUN
4167 IF (IWATER(I).EQ.1) VX0PSIM(I) = VAPSIM(I)
4168 9000 CONTINUE
4169 ENDIF
4170 IF ( ITYPE .GE. 3 ) THEN
4171 DO 9002 I = 1,IRUN
4172 VX0PSIM(I) = VX0(I) * VAPSIM(I)
4173 9002 CONTINUE
4174 ENDIF
4175 IF ( ITYPE .NE. 2 ) THEN
13fa5cb63c Jean*4176
1662f365b2 Andr*4177 DO 9004 I = 1,IRUN
4178 VDZ0(I) = 0.
4179 VG(I) = 0.
4180 VG0(I) = 0.
4181 VR1MG0(I) = 1.
4182 9004 CONTINUE
13fa5cb63c Jean*4183
1662f365b2 Andr*4184 IF ( LWATER ) THEN
4185 CALL ZCSUB ( VUSTAR,VDZSEA,IWATER,.TRUE.,IRUN,VZ2)
13fa5cb63c Jean*4186
1662f365b2 Andr*4187 DO 9006 I = 1,IRUN
4188 IF ( IWATER(I).EQ.1) THEN
4189 VAZ0(I) = 1. / VZ1(I)
4190 VG(I) = VDZSEA(I) * VAZ0(I)
4191 VG0(I) = VX0PSIM(I) * VG(I)
4192 VR1MG0(I) = 1. / ( 1. - VG0(I) )
4193 VDZ0(I) = ( VZ2(I) - VZ1(I) ) * VR1MG0(I)
4194 ENDIF
4195 9006 CONTINUE
4196 ENDIF
4197 ENDIF
13fa5cb63c Jean*4198
1662f365b2 Andr*4199 IF ( LWATER .AND. (ITYPE.GE.3) ) THEN
4200 DO 9008 I = 1,IRUN
4201 IF (IWATER(I).EQ.1) VDZ0(I) = VDZ0(I) * VAZ0(I)
4202 9008 CONTINUE
4203 ENDIF
13fa5cb63c Jean*4204
1662f365b2 Andr*4205
13fa5cb63c Jean*4206
1662f365b2 Andr*4207 IF (ITYPE.GE.3) THEN
4208 DO 9010 I = 1,IRUN
4209 VXNUM1(I) = 0.
4210 IF (VRIB1(I).EQ.0.) THEN
4211 INTRIB(I) = 1
4212 ELSE
4213 INTRIB(I) = 0
4214 ENDIF
4215 IF ( INTRIB(I).EQ.0 ) VXNUM1(I) = 1. / VRIB1(I)
4216 VPSIGB2(I) = 0.
4217 if(vpsig(i).gt.0.)VPSIGB2(I) =
4218 1 0.5 * ( vpsig(i)*vpsig(i) + b2uhs ) / vpsig(i)
4219 VDX(I) = VX(I) - VX0(I)
4220 VDXPSIM(I) = VDX(I) * VAPSIM(I)
4221 VDY(I) = VY(I) - VY0(I)
4222 VXNUM3(I) = - VPSIGB2(I)
13fa5cb63c Jean*4223
1662f365b2 Andr*4224 IF ( LWATER ) THEN
4225
4226 IF (IWATER(I).EQ.1) THEN
4227 VDXPSIM(I) = VDXPSIM(I) * VR1MG0(I)
4228 VXNUM3(I) = VXNUM3(I) + VG(I) * ( VY0(I) - VPSIGB2(I) )
4229 VXNUM2(I) = VY0(I) - VPSIGB2(I) - VX0PSIM(I) * VPSIGB2(I)
4230 VXNUM2(I) = (VXNUM2(I) * VAPSIHG(I)) - 2. * VX0PSIM(I)
4231 VXNUM2(I) = VXNUM2(I) * VDZ0(I)
4232 ENDIF
4233 ENDIF
13fa5cb63c Jean*4234
1662f365b2 Andr*4235 VDEN(I) = VDY(I) + VDXPSIM(I) * VXNUM3(I)
4236 VDEN(I) = ( 1. + VDEN(I) * VAPSIHG(I) ) - 2. * VDXPSIM(I)
4237 9010 CONTINUE
4238 ENDIF
13fa5cb63c Jean*4239
1662f365b2 Andr*4240 IF (ITYPE.EQ.5) THEN
4241 DO 9012 I = 1,IRUN
4242 VAWS1(I) = VR1MG0(I) / VWS1(I)
4243 VXNUM3(I) = VXNUM3(I) * VAPSIHG(I)
13fa5cb63c Jean*4244
1662f365b2 Andr*4245 IF ( LWATER ) THEN
4246
4247 IF(IWATER(I).EQ.1) THEN
4248 VXNUM3(I) = VXNUM3(I) - 2. * VG0(I)
4249 VXNUM3(I) = VAWS1(I) * VXNUM3(I)
4250 ENDIF
4251 ENDIF
4252 9012 CONTINUE
4253 ENDIF
13fa5cb63c Jean*4254
1662f365b2 Andr*4255
13fa5cb63c Jean*4256
1662f365b2 Andr*4257 IF (ITYPE.GE.2) THEN
4258 DO 9014 I = 1,IRUN
4259 VXNUM(I) = VRIB2(I) - VRIB1(I)
4260 IF( (VX0(I).GT.XX0MAX).AND.(VXNUM(I).GE.0.) )VXNUM(I) = 0.
4261 VXNUM(I) = VXNUM1(I) * VXNUM(I)
4262 9014 CONTINUE
4263 ENDIF
13fa5cb63c Jean*4264
1662f365b2 Andr*4265 IF ( ITYPE.EQ.2 )THEN
4266 DO 9016 I = 1,IRUN
4267 VDZETA1(I) = VDZETA(I)
4268 VXNUM(I) = VXNUM(I) + VXNUM3(I) * ( VWS2(I) - VWS1(I) )
4269 9016 CONTINUE
4270 ENDIF
13fa5cb63c Jean*4271
1662f365b2 Andr*4272 IF (ITYPE.GE.3) THEN
4273 DO 9018 I = 1,IRUN
4274 VDZETA1(I) = VXNUM(I)
4275 IF(LWATER.AND.(IWATER(I).EQ.1)) VXNUM(I) = VXNUM(I) + VXNUM2(I)
4276 IF ( VDEN(I) .LT.0.1 ) VDEN(I) = 0.1
4277 9018 CONTINUE
4278 ENDIF
13fa5cb63c Jean*4279
1662f365b2 Andr*4280 IF (ITYPE.GE.2) THEN
4281 DO 9020 I = 1,IRUN
4282 VDZETA(I) = VXNUM(I) / VDEN(I)
4283 9020 CONTINUE
4284 ENDIF
4285 IF (ITYPE.GE.3) THEN
4286 DO 9022 I = 1,IRUN
4287 IF ( (VRIB2(I).EQ.0.) .OR. (VDZETA(I).LE.-1.) ) THEN
4288 VINT1(I) = 1
4289 ELSE
4290 VINT1(I) = 0
4291 ENDIF
4292 IF ( VINT1(I).EQ.1 ) VDZETA(I) = VDZETA1(I)
4293 9022 CONTINUE
4294 ENDIF
4295 IF (ITYPE.EQ.2) THEN
4296 DO 9024 I = 1,IRUN
4297 VDZETA2(I) = VDZETA(I) + VDZETA1(I)
4298 IF ( (VRIB2(I).EQ.0.) .OR. (VDZETA2(I).LE.-1.) ) THEN
4299 VINT1(I) = 1
4300 ELSE
4301 VINT1(I) = 0
4302 ENDIF
4303 IF(VINT1(I).EQ.1)VDZETA(I)=VXNUM1(I)*VRIB2(I) - 1. - VDZETA1(I)
4304 9024 CONTINUE
4305 ENDIF
4306
4307
13fa5cb63c Jean*4308
1662f365b2 Andr*4309 IF ( LWATER .AND. (ITYPE.GE.3) )THEN
4310 DO 9026 I = 1,IRUN
4311 IF( IWATER(I).EQ.1 ) THEN
4312 VZCOEF2(I) = VG(I) * VDXPSIM(I)
4313 VDZ0(I) = VDZ0(I) - VZCOEF2(I) * VDZETA(I)
4314 ENDIF
4315 9026 CONTINUE
4316 ENDIF
13fa5cb63c Jean*4317
1662f365b2 Andr*4318 IF ( LWATER .AND. (ITYPE.EQ.5) ) THEN
4319 DO 9028 I = 1,IRUN
4320 IF(IWATER(I).EQ.1) VZCOEF1(I) = VG(I) * VAWS1(I)
4321 9028 CONTINUE
4322 ENDIF
13fa5cb63c Jean*4323
1662f365b2 Andr*4324 IF ( LWATER .AND. (ITYPE.EQ.2) ) THEN
4325 DO 9030 I = 1,IRUN
4326 IF (IWATER(I).EQ.1) VDZ0(I) =
4327 1 VZCOEF1(I) * ( VWS2(I) - VWS1(I) ) - VZCOEF2(I) * VDZETA(I)
4328 9030 CONTINUE
4329 ENDIF
13fa5cb63c Jean*4330
1662f365b2 Andr*4331
13fa5cb63c Jean*4332
1662f365b2 Andr*4333 IF ( (ITYPE.EQ.1) .AND. LWATER ) THEN
4334 DO 9032 I = 1,IRUN
4335 IF (IWATER(I).EQ.1) THEN
4336 VDPSIM(I) = - VDZ0(I) * VAZ0(I)
4337 VDPSIH(I) = VDPSIM(I)
4338 ENDIF
4339 9032 CONTINUE
4340 ENDIF
13fa5cb63c Jean*4341
1662f365b2 Andr*4342 IF (ITYPE.GE.3) THEN
4343 DO 9034 I = 1,IRUN
4344 VDPSIM(I) = VDX(I) * VDZETA(I)
4345 VDPSIH(I) = VDY(I) * VDZETA(I)
4346 IF ( LWATER ) THEN
4347 IF (IWATER(I).EQ.1 ) THEN
4348 VDPSIM(I) = VDPSIM(I) - VX0(I) * VDZ0(I)
4349 VDPSIH(I) = VDPSIH(I) - VY0(I) * VDZ0(I)
4350 ENDIF
4351 ENDIF
4352 9034 CONTINUE
4353 ENDIF
13fa5cb63c Jean*4354
1662f365b2 Andr*4355
13fa5cb63c Jean*4356
1662f365b2 Andr*4357 IF (ITYPE.GE.4) THEN
4358 DO 9036 I = 1,IRUN
4359 VDPSIMC(I) = -0.9 - VDPSIM(I) * VAPSIM(I)
4360 VDPSIHC(I) = -0.9 * VPSIH(I) - VDPSIH(I)
4361 IF ( VDPSIMC(I).GT.0. ) THEN
4362 VINT1(I) = 1
4363 ELSE
4364 VINT1(I) = 0
4365 ENDIF
4366 IF ( VDPSIHC(I).GT.0. ) THEN
4367 VINT2(I) = 1
4368 ELSE
4369 VINT2(I) = 0
4370 ENDIF
4371 VDZETA1(I) = 0.
4372 IF(VINT1(I).EQ.1) VDZETA1(I) = VDPSIMC(I) / VDXPSIM(I)
4373 IF((VINT1(I).EQ.1).OR.(VINT2(I).EQ.1)) VTEMPLIN(I) =
4374 1 VDY(I) + VY0(I) * VG(I) * VDXPSIM(I)
4375 IF (VINT2(I).EQ.1) then
4376 VDZETA2(I) = VDPSIHC(I) / VTEMPLIN(I)
4377 IF ( VDZETA2(I).LT.VDZETA1(I) ) VDZETA1(I) = VDZETA2(I)
4378 endif
4379 IF((VINT1(I).EQ.1).OR.(VINT2(I).EQ.1)) THEN
4380 VDZETA(I) = VDZETA1(I) + VDZETA(I)
4381 VDPSIM(I) = VDPSIM(I) + VDX(I) * VR1MG0(I) * VDZETA1(I)
4382 VDPSIH(I) = VDPSIH(I) + VTEMPLIN(I) * VDZETA1(I)
4383 IF ( IWATER(I).EQ.1 )
4384 1 VDZ0(I) = VDZ0(I) - VG(I) * VDXPSIM(I) * VDZETA1(I)
4385 ENDIF
4386 9036 CONTINUE
4387 ENDIF
13fa5cb63c Jean*4388
1662f365b2 Andr*4389 RETURN
4390 END
13fa5cb63c Jean*4391
4392
4393
1662f365b2 Andr*4394 SUBROUTINE ZCSUB (VUSTAR,VDZSEA,IWATER,LDZSEA,IRUN,VZSEA)
4395
4396
4397
4398
4399
4400
4401
4402
4403
4404
4405
4406
4407
4408
4409
4410
4411
4412
4413
4414
5317312052 Jean*4415 implicit none
453e05dab3 Andr*4416
4417
4418 integer irun
a456aa407c Andr*4419 _RL VZSEA(IRUN),VUSTAR(IRUN),VDZSEA(IRUN)
453e05dab3 Andr*4420 integer IWATER(IRUN)
4421 LOGICAL LDZSEA
4422
4423
a456aa407c Andr*4424 _RL USTMX1,USTMX2,USTMX3
1662f365b2 Andr*4425 PARAMETER ( USTMX1 = 1.14973 )
4426 PARAMETER ( USTMX2 = 0.381844 )
4427 PARAMETER ( USTMX3 = 0.0632456)
4428
a456aa407c Andr*4429 _RL AA(IRUN,5),TEMP(IRUN)
453e05dab3 Andr*4430 integer INT2(IRUN),INT3(IRUN),INT4(IRUN)
4431 integer i,k
4432
a456aa407c Andr*4433 _RL AA1(5),AA2(5),AA3(5),AA4(5)
1662f365b2 Andr*4434 DATA AA1/.2030325E-5,0.0,0.0,0.0,0.0/
4435 DATA AA2/-0.402451E-08,0.239597E-04,0.117484E-03,0.191918E-03,
4436 1 0.395649E-04/
4437 DATA AA3/-0.237910E-04,0.228221E-03,-0.860810E-03,0.176543E-02,
4438 1 0.784260E-04/
4439 DATA AA4/-0.343228E-04,0.552305E-03,-0.167541E-02,0.250208E-02,
4440 1 -0.153259E-03/
13fa5cb63c Jean*4441
1662f365b2 Andr*4442
4443
4444
13fa5cb63c Jean*4445
1662f365b2 Andr*4446 DO 9000 I = 1,IRUN
4447 IF(VUSTAR(I) .LT. 1.e-6)THEN
4448 INT3(I) = 1
4449 ELSE
4450 INT3(I) = 0
4451 ENDIF
4452 9000 CONTINUE
4453 DO 9002 I = 1,IRUN
4454 IF(INT3(I).EQ.1) VUSTAR(I) = 1.e-6
4455 9002 CONTINUE
13fa5cb63c Jean*4456
1662f365b2 Andr*4457
4458
4459
13fa5cb63c Jean*4460
1662f365b2 Andr*4461 DO 9004 I = 1,IRUN
4462 IF( (VUSTAR(I) .GT. USTMX1) .AND. (IWATER(I).EQ.1) ) THEN
4463 INT4(I) = 1
4464 ELSE
4465 INT4(I) = 0
4466 ENDIF
4467 9004 CONTINUE
4468 DO 9006 I = 1,IRUN
4469 IF(VUSTAR(I) .GT. USTMX2) THEN
4470 INT3(I) = 1
4471 ELSE
4472 INT3(I) = 0
4473 ENDIF
4474 9006 CONTINUE
4475 DO 9008 I = 1,IRUN
4476 IF(VUSTAR(I) .GE. USTMX3) THEN
4477 INT2(I) = 1
4478 ELSE
4479 INT2(I) = 0
4480 ENDIF
4481 9008 CONTINUE
13fa5cb63c Jean*4482
1662f365b2 Andr*4483 DO 100 K=1,5
4484 DO 9010 I = 1,IRUN
4485 AA(I,K) = AA1(K)
4486 IF( INT2(I).EQ.1 ) AA(I,K) = AA2(K)
4487 IF( INT3(I).EQ.1 ) AA(I,K) = AA3(K)
4488 IF( INT4(I).EQ.1 ) AA(I,K) = AA4(K)
4489 9010 CONTINUE
4490 100 CONTINUE
13fa5cb63c Jean*4491
1662f365b2 Andr*4492
4493
4494
13fa5cb63c Jean*4495
1662f365b2 Andr*4496 DO 9012 I = 1,IRUN
4497 VDZSEA(I) = ( AA(I,4) + AA(I,5) * VUSTAR(I) ) * VUSTAR(I)
4498 VZSEA(I) = AA(I,2) + ( AA(I,3) + VDZSEA(I) ) * VUSTAR(I)
4499 TEMP(I) = AA(I,1) / VUSTAR(I)
4500 VZSEA(I) = VZSEA(I) + TEMP(I)
4501 9012 CONTINUE
13fa5cb63c Jean*4502
1662f365b2 Andr*4503
4504
4505
13fa5cb63c Jean*4506
1662f365b2 Andr*4507 IF( LDZSEA ) THEN
4508 DO 9014 I = 1,IRUN
4509 VDZSEA(I) = 3. * VDZSEA(I) -(AA(I,4)*VUSTAR(I) - AA(I,3))
4510 VDZSEA(I) = VDZSEA(I) * VUSTAR(I) - TEMP(I)
4511 9014 CONTINUE
4512 ENDIF
13fa5cb63c Jean*4513
1662f365b2 Andr*4514 RETURN
4515 END
4516
13fa5cb63c Jean*4517
4518
4519 SUBROUTINE SEAICE ( nocean, timstp, hice,
4520 & eturb, dedtc,
4521 & hsturb, dhsdtc,
4522 & qice, dqice,
4523 & swnet, lwnet, dst4,
4524 & pke, seaic, tc, qa )
1662f365b2 Andr*4525 implicit none
4526 integer nocean
a456aa407c Andr*4527 _RL timstp
4528 _RL eturb(nocean),dedtc(nocean),hsturb(nocean),dhsdtc(nocean)
4529 _RL swnet(nocean),lwnet(nocean), dst4(nocean)
4530 _RL qice(nocean),dqice(nocean)
4531 _RL pke(nocean), tc(nocean), qa(nocean)
4532 _RL seaic(nocean)
1662f365b2 Andr*4533
4534
a456aa407c Andr*4535 _RL rhoC
453e05dab3 Andr*4536 parameter (rhoC = 1.93e6)
1662f365b2 Andr*4537
a456aa407c Andr*4538 _RL faceps,getcon,latent,codt,deltg,hice
1662f365b2 Andr*4539 integer i
4540
4541 faceps = getcon('EPSFAC')
4542 latent = getcon('HEATI') * getcon('CALTOJ')
4543
5317312052 Jean*4544 codt = rhoC * (hice/100) / timstp
1662f365b2 Andr*4545
5317312052 Jean*4546
4547
1662f365b2 Andr*4548 do i =1,nocean
4549 if( seaic(i).gt.0.0 ) then
4550 deltg = ( swnet(i)-lwnet(i)-latent*eturb(i)-hsturb(i)+qice(i) )
13fa5cb63c Jean*4551 & / ( codt+dst4(i)+latent*dedtc(i)+dhsdtc(i)-dqice(i) )
1662f365b2 Andr*4552 qa(i) = qa(i) + (faceps*qa(i)/(tc(i)*tc(i)))*deltg
5317312052 Jean*4553 tc(i) = tc(i) + deltg
1662f365b2 Andr*4554 endif
4555 enddo
4556
4557 return
4558 end