45subroutine rtsig(lusig,head,k1,k2,kgds,ijo,levs,ntrac,jcap,lnt2,me, &
46 h,p,px,py,t,u,v,d,trc,iret)
48 use sigio_module,
only : sigio_intkind, sigio_head
49 use sigio_r_module,
only : sigio_dati, sigio_rrdati
50 use physcons_post,
only : con_omega, con_fvirt
53 integer(sigio_intkind),
intent(in) :: lusig
54 type(sigio_head),
intent(in) :: head
55 integer,
intent(in) :: k1,k2,kgds(200),ijo,levs,ntrac,jcap,lnt2,me
56 real,
dimension(ijo),
intent(out) :: h,p,px,py
57 real,
dimension(ijo,k1:k2),
intent(out):: t,u,v,d
58 real,
dimension(ijo,k1:k2,ntrac),
intent(out),
target :: trc
59 integer,
intent(out) :: iret
61 integer idrt,io,jo,iidea
63 integer(sigio_intkind):: irets
70 real,
target :: trisca(lnt2,k1:k2+1), triscb(lnt2,k1:k2)
71 real,
dimension(:),
allocatable :: cpi
72 real,
dimension(:,:),
allocatable :: wrk
73 integer io3,ict,jct,n,i,k,jc,nt
80 if(kgds(1) == 0 .and. kgds(4) < 90000) idrt = 256
87 print*,
'Debug rtsig= ',lusig,k1,k2,ijo,kgds(1:20)
98 dati%f => trisca(:,k1)
99 call sigio_rrdati(lusig,head,dati,irets)
100 if(irets /= 0)
return
108 dati%f => trisca(:,k1+1)
109 call sigio_rrdati(lusig,head,dati,irets)
110 if(irets /= 0)
return
116 call sptezm(0,jcap,idrt,io,jo,2,trisca(1,k1),wrk,1)
117 if( mod(idvm,10) < 2)
then
121 p(i) = 1.e3*exp(wrk(i,2))
124 elseif(mod(idvm,10) == 2)
then
128 p(i) = 1000.*wrk(i,2)
132 if (
allocated(wrk))
deallocate(wrk)
134 call sptezd(0,jcap,idrt,io,jo,trisca(1,k1+1),pmean,px,py,1)
145 write(*,*)
' retriving T for k=',k,
' k1=',k1,
' k2=',k2
147 dati%f => trisca(:,k)
149 call sigio_rrdati(lusig,head,dati,iret)
151 call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),t(1,k1),1)
154 dati%i = levs + 2 + (k-1) * 2 + 1
155 dati%f => trisca(:,k)
156 call sigio_rrdati(lusig,head,dati,irets)
157 if(irets /= 0)
return
158 dati%i = levs + 2 + (k-1) * 2 + 2
159 dati%f => triscb(:,k)
160 call sigio_rrdati(lusig,head,dati,irets)
161 if(irets /= 0)
return
163 call sptezmv(0,jcap,idrt,io,jo,klen,trisca(1,k1),triscb(1,k1), &
165 call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),d(1,k1),1)
168 write(*,*)
' retriving d/z for k=',k,
' k1=',k1,
' k2=',k2
171 write(*,*)
' start get tracer'
174 dati%i = levs * (2+nt) + 2 + k
175 dati%f => trisca(:,k)
176 call sigio_rrdati(lusig,head,dati,irets)
178 call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),trc(1,k1,nt),1)
179 write(*,*)
' retriving d/z for nt=',nt,
'ntrac=',ntrac,
'k=',k,
' k1=',k1,
' k2=',k2
182 write(*,*)
' end get tracer,idvm=',idvm,
'ijo=',ijo,
'ntrac=',ntrac
185 if (mod(idvm/10,10) == 3)
then
186 allocate(cpi(0:ntrac))
188 cpi(0:ntrac) = head%cpi(1:ntrac+1)
196 if( cpi(n) /= 0.0 )
then
197 xcp = xcp + cpi(n)*trc(i,k,n)
198 sumq = sumq + trc(i,k,n)
201 xcp = (1.-sumq)*cpi(0) + xcp
202 t(i,k) = t(i,k) / xcp
205 if (
allocated(cpi))
deallocate(cpi)
210 t(i,k) = t(i,k) / (1+con_fvirt*trc(i,k,1))
327 subroutine modstuff2(im,ix,km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,&
329 use sigio_module,
only : sigio_modprd
331 integer,
intent(in) :: im,ix,km,idvc,idsl,nvcoord,me
332 real,
intent(in) :: vcoord(km+1,nvcoord)
333 real,
dimension(ix),
intent(in) :: ps,psx,psy
334 real,
dimension(ix,km),
intent(in) :: u,v,d
335 real*8,
dimension(ix,km+1),
intent(out) :: pi
336 real*8,
dimension(ix,km),
intent(out) :: pm
337 real,
dimension(ix,km),
intent(out) :: om
345 real*8 ps8(ix),pm8(ix,km),pd8(ix,km),vcoord8(km+1,nvcoord)
346 real*8 dpmdps(ix,km),dpddps(ix,km),dpidps(ix,km+1),pi8(ix,km+1)
347 real vgradpps,pd(im,km),os(im)
349 integer i,k,iret,logk
354 call sigio_modprd(im,ix,km,nvcoord,idvc,idsl,vcoord8,iret, &
355 ps=ps8,pd=pd8,dpddps=dpddps,pm=pm8,dpmdps=dpmdps)
374 pi8(i,k+1) = pi8(i,k) - pd8(i,k)
375 dpidps(i,k+1) = dpidps(i,k) - dpddps(i,k)
379 pi(i,k+1) = pi8(i,k+1)
387 vgradpps = (u(i,k)*psx(i) + v(i,k)*psy(i)) * ps(i)
389 os(i) = os(i) - vgradpps*(dpmdps(i,k)-dpidps(i,k+1)) &
390 - d(i,k)*(pm(i,k)-pi(i,k+1))
392 om(i,k) = os(i) + vgradpps*dpmdps(i,k)
394 os(i) = os(i) - vgradpps*(dpidps(i,k)-dpmdps(i,k)) &
395 - d(i,k)*(pi(i,k)-pm(i,k))
448 subroutine trssc(jcap,nc,km,ntrac,idvc,idvm,idsl,nvcoord,vcoord, &
449 cpi,idrt,lonb,latb,ijl,ijn,j1,j2,jc,chgq0, &
450 szs,sps,st,sd,sz,sq,gfszs,gfsps,gfsp,gfsdp, &
451 gfst,gfsu,gfsv,gfsq,gfsw)
453 integer,
intent(in)::jcap,nc,km,ntrac,idvc,idvm,idsl,nvcoord,idrt,lonb,latb
454 integer,
intent(in)::ijl,ijn,j1,j2,jc,chgq0
455 real,
intent(in):: szs(nc),sps(nc),st(nc,km),sd(nc,km),sz(nc,km),sq(nc,km*ntrac)
456 real,
intent(in):: cpi(0:ntrac)
457 real*8,
intent(in):: vcoord(km+1,nvcoord)
458 real,
dimension(ijn),
intent(inout):: gfszs,gfsps
459 real,
dimension(ijn,km),
intent(inout):: gfsp,gfsdp,gfst,gfsu,gfsv,gfsw
460 real,
dimension(ijn,km*ntrac),
intent(inout):: gfsq
461 real zs(ijl),ps(ijl),t(ijl,km),u(ijl,km),v(ijl,km),q(ijl,km*ntrac)
462 real wi(ijl,km),pi(ijl,km),dpo(ijl,km)
464 integer thermodyn_id,jn,js,is,in
465 integer jj,jjm,k,n,j,i,ij,lonb2
468 if(j1==732)print*,
'sample input to trssc= ',jcap,nc,km,ntrac, &
469 idvc,idvm,idsl,nvcoord, &
470 idrt,lonb,latb,ijl,ijn,j1,j2,jc,chgq0
475 call sptran(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijl, &
476 j1,j2,jc,szs,zs(in),zs(is),1)
477 call sptran(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijl, &
478 j1,j2,jc,sps,ps(in),ps(is),1)
479 call sptran(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijl, &
480 j1,j2,jc,st,t(in,1),t(is,1),1)
481 call sptranv(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijl, &
482 j1,j2,jc,sd,sz,u(in,1),u(is,1),v(in,1),v(is,1),1)
483 call sptran(0,jcap,idrt,lonb,latb,km*ntrac,1,1,lonb2,lonb2,nc,ijl, &
484 j1,j2,jc,sq,q(in,1),q(is,1),1)
488 if(t(i,k)>400. .or. t(i,k)<100.)print*,
'bad T from sptran',i,k,t(i,k)
489 if(q(i,k)>1.)print*,
'bad Q from sptran',i,k,q(i,k)
490 if(q(i,2*k)>1.)print*,
'bad Q from sptran',i,k,q(i,2*k)
491 if(q(i,3*k)>1.)print*,
'bad Q from sptran',i,k,q(i,3*k)
495 select case(mod(idvm,10))
498 ps(i)=1.e3*exp(ps(i))
506 ps(i)=1.e3*exp(ps(i))
510 thermodyn_id=mod(idvm/10,10)
511 if (thermodyn_id == 3)
then
514 t(i,k) = t(i,k)/cpi(0)
520 call getomega(jcap,nc,km,idvc,idvm,idrt,idsl, &
521 nvcoord,vcoord,lonb,latb,ijl,j1,j2,1,sd,sps, &
526 if(t(i,k)>400. .or. t(i,k)<100.)print*,
'bad T after getomega',i,k,t(i,k)
527 if(q(i,k)>1. )print*,
'bad Q after getomega',i,k,q(i,k)
528 if(q(i,2*k)>1. )print*,
'bad Q after getomega',i,2*k,q(i,2*k)
532 if(thermodyn_id /= 2)
then
534 if (thermodyn_id == 3)
then
540 if( cpi(n) /= 0.0 )
then
541 xcp = xcp + cpi(n)*q(i,k+(n-1)*km)
542 sumq = sumq + q(i,k+(n-1)*km)
545 t(i,k) = t(i,k)/((1.-sumq)*cpi(0)+xcp)
550 tvcon=(461.50/287.05-1.)
551 t(:,:) = t(:,:)/(1.+tvcon*q(:,1:km))
557 if(t(i,k)>400. .or. t(i,k)<100.)print*,
'bad T after Tv to T',i,k,t(i,k)
558 if(q(i,k)>1.)print*,
'bad Q after Tv to T',i,k,q(i,k)
559 if(q(i,2*k)>1. )print*,
'bad Q after Tv to T',i,k,q(i,2*k)
565 if (chgq0 == 1) q = max(q, 0.0)
576 gfszs(i+jn) = zs(i+jjm)
577 gfsps(i+jn) = ps(i+jjm)
580 gfszs(i+js) = zs(i+jj)
581 gfsps(i+js) = ps(i+jj)
585 gfsdp(i+jn,k) = dpo(i+jjm,k)
586 gfsp(i+jn,k) = pi(i+jjm,k)
587 gfst(i+jn,k) = t(i+jjm,k)
588 gfsu(i+jn,k) = u(i+jjm,k)
589 gfsv(i+jn,k) = v(i+jjm,k)
590 gfsw(i+jn,k) = wi(i+jjm,k)
593 gfsdp(i+js,k) = dpo(i+jj,k)
594 gfsp(i+js,k) = pi(i+jj,k)
595 gfst(i+js,k) = t(i+jj,k)
596 gfsu(i+js,k) = u(i+jj,k)
597 gfsv(i+js,k) = v(i+jj,k)
598 gfsw(i+js,k) = wi(i+jj,k)
603 gfsq(i+jn,k) = q(i+jjm,k)
606 gfsq(i+js,k) = q(i+jj,k)
640 subroutine getomega(jcap,nc,km,idvc,idvm,idrt,idsl,nvcoord,vcoord, &
641 lonb,latb,ijn,j1,j2,jc,sd,sps,psi,ti,ui,vi,wi,pm,pd)
643 use sigio_module,
only : sigio_modprd
646 integer,
intent(in):: jcap,nc,km,idvc,idvm,idrt,idsl,nvcoord
647 integer,
intent(in):: lonb,latb,j1,j2,jc,ijn
648 real*8,
intent(in):: vcoord(km+1,nvcoord)
649 real,
intent(in):: sd(nc,km),sps(nc)
650 real,
intent(in):: psi(ijn),ti(ijn,km),ui(ijn,km),vi(ijn,km)
651 real,
intent(out):: wi(ijn,km),pm(ijn,km),pd(ijn,km)
654 real*8 psi8(ijn),ti8(ijn,km),pm8(ijn,km),pd8(ijn,km)
655 real*8 dpmdps(ijn,km),dpddps(ijn,km),dpidps(ijn,km+1),vgradp,psmean
656 real di(ijn,km),psx(ijn),psy(ijn)
657 integer k,i,ij,lonb2,iret,is,in
663 call sptrand(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijn, &
664 j1,j2,jc,sps,psmean,psx(in),psx(is),psy(in),psy(is),1)
666 call sptran(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijn, &
667 j1,j2,jc,sd,di(in,1),di(is,1),1)
671 call sigio_modprd(ijn,ijn,km,nvcoord,idvc,idsl,vcoord,iret, &
672 ps=psi8,t=ti8,pm=pm8,pd=pd8,dpmdps=dpmdps,dpddps=dpddps)
676 select case(mod(idvm,10))
681 psx(i)=psx(i)/(psi(i)*1.0e-3)
682 psy(i)=psy(i)/(psi(i)*1.0e-3)
698 pi(i,k+1)=pi(i,k)-pd(i,k)
699 dpidps(i,k+1)=dpidps(i,k)-dpddps(i,k)
705 vgradp=ui(i,k)*psx(i)+vi(i,k)*psy(i)
706 os=os-vgradp*psi(i)*(dpmdps(i,k)-dpidps(i,k+1))- &
707 di(i,k)*(pm(i,k)-pi(i,k+1))
708 wi(i,k)=vgradp*psi(i)*dpmdps(i,k)+os
709 os=os-vgradp*psi(i)*(dpidps(i,k)-dpmdps(i,k))- &
710 di(i,k)*(pi(i,k)-pm(i,k))
subroutine modstuff(km, idvc, idsl, nvcoord, vcoord, ps, psx, psy, d, u, v, pi, pm, om)
Computes model coordinate dependent functions.
subroutine rtsig(lusig, head, k1, k2, kgds, ijo, levs, ntrac, jcap, lnt2, me, h, p, px, py, t, u, v, d, trc, iret)
Reads and transforms a sigma file.
subroutine trssc(jcap, nc, km, ntrac, idvc, idvm, idsl, nvcoord, vcoord, cpi, idrt, lonb, latb, ijl, ijn, j1, j2, jc, chgq0, szs, sps, st, sd, sz, sq, gfszs, gfsps, gfsp, gfsdp, gfst, gfsu, gfsv, gfsq, gfsw)
Transforms sigma spectral fields to grid.
subroutine getomega(jcap, nc, km, idvc, idvm, idrt, idsl, nvcoord, vcoord, lonb, latb, ijn, j1, j2, jc, sd, sps, psi, ti, ui, vi, wi, pm, pd)
Computes omega.
subroutine modstuff2(im, ix, km, idvc, idsl, nvcoord, vcoord, ps, psx, psy, d, u, v, pi, pm, om, me)
Computes model coordinate dependent functions.