UPP  11.0.0
 All Data Structures Files Functions Variables Pages
GFSPOSTSIG.F
Go to the documentation of this file.
1 
50 
51 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94 subroutine rtsig(lusig,head,k1,k2,kgds,ijo,levs,ntrac,jcap,lnt2,me, &
95  h,p,px,py,t,u,v,d,trc,iret)
96 
97  use sigio_module, only : sigio_intkind, sigio_head
98  use sigio_r_module, only : sigio_dati, sigio_rrdati
99  use physcons_post, only : con_omega, con_fvirt
100  use omp_lib
101  implicit none
102  integer(sigio_intkind),intent(in) :: lusig
103  type(sigio_head), intent(in) :: head
104  integer,intent(in) :: k1,k2,kgds(200),ijo,levs,ntrac,jcap,lnt2,me
105  real,dimension(ijo), intent(out) :: h,p,px,py
106  real,dimension(ijo,k1:k2),intent(out):: t,u,v,d
107  real,dimension(ijo,k1:k2,ntrac),intent(out),target :: trc
108  integer,intent(out) :: iret
109 !
110  integer idrt,io,jo,iidea
111 ! integer idrt,io,jo,mo3,mct,iidea,midea
112  integer(sigio_intkind):: irets
113 ! type(sigio_datm):: datm
114  type(sigio_dati) dati
115 ! type griddata
116 ! real,dimension(:,:),pointer :: datm
117 ! endtype griddata
118 ! type(griddata),dimension(:),pointer :: datatrc
119  real, target :: trisca(lnt2,k1:k2+1), triscb(lnt2,k1:k2)
120  real,dimension(:), allocatable :: cpi
121  real,dimension(:,:),allocatable :: wrk
122  integer io3,ict,jct,n,i,k,jc,nt
123  integer idvm, klen
124  real pmean,sumq,xcp
125 ! integer, parameter :: latch=20
126 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127 ! Determine output grid
128  idrt = kgds(1)
129  if(kgds(1) == 0 .and. kgds(4) < 90000) idrt = 256
130  io = kgds(2)
131  jo = kgds(3)
132 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133 ! Read and transform surface fields
134  iret = 1
135  if (me == 0) then
136  print*,'Debug rtsig= ',lusig,k1,k2,ijo,kgds(1:20)
137  endif
138 
139  idvm = head%idvm
140 ! jc = omp_get_num_threads()
141 ! write(*,*)' in RTSIG lnt2=',lnt2,' threads=',jc,' latch=',latch, &
142 ! ' jcap=',jcap,' io=',io,' jo=',jo,' ijo=',ijo
143 !
144  if (k2 < k1) return
145 
146  dati%i = 1 ! hs
147  dati%f => trisca(:,k1)
148  call sigio_rrdati(lusig,head,dati,irets)
149  if(irets /= 0) return
150 
151 ! call sptez(0,jcap,idrt,io,jo,trisca(1,k1),h,1)
152 ! call sptez(0,jcap,idrt,io,jo,dats%hs,h,1)
153 ! call sptez(0,jcap,idrt,io,jo,dats%ps,p,1)
154 ! call sptezj(jcap,lnt2,1,idrt,io,jo,jc,trisca,h,latch,1)
155 !
156  dati%i = 2 ! Surface pressure
157  dati%f => trisca(:,k1+1)
158  call sigio_rrdati(lusig,head,dati,irets)
159  if(irets /= 0) return
160 !
161 ! call sptez(0,jcap,idrt,io,jo,trisca(1,k1),p,1)
162 ! call sptezj(jcap,lnt2,1,idrt,io,jo,jc,trisca,p,latch,1)
163 !--
164  allocate(wrk(ijo,2))
165  call sptezm(0,jcap,idrt,io,jo,2,trisca(1,k1),wrk,1)
166  if( mod(idvm,10) < 2) then
167 !$omp parallel do private(i)
168  do i=1,ijo
169  h(i) = wrk(i,1)
170  p(i) = 1.e3*exp(wrk(i,2))
171 ! p(i) = 1.e3*exp(p(i))
172  enddo
173  elseif(mod(idvm,10) == 2) then
174 !$omp parallel do private(i)
175  do i=1,ijo
176  h(i) = wrk(i,1)
177  p(i) = 1000.*wrk(i,2)
178 ! p(i) = 1000.*p(i)
179  enddo
180  endif
181  if (allocated(wrk)) deallocate(wrk)
182 
183  call sptezd(0,jcap,idrt,io,jo,trisca(1,k1+1),pmean,px,py,1)
184  iret = 0
185 
186 ! if (k2 < k1) return
187 
188 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189 ! Read and transform fields on levels k1 through k2
190  iret = 2
191  if (k2 >= k1) then
192  klen = k2-k1+1
193  do k=k1,k2
194  write(*,*)' retriving T for k=',k,' k1=',k1,' k2=',k2
195  dati%i = k + 2 ! Virtual Temperature or CpT
196  dati%f => trisca(:,k)
197 
198  call sigio_rrdati(lusig,head,dati,iret)
199  enddo
200  call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),t(1,k1),1)
201 ! call sptezm(0,jcap,idrt,io,jo,klen,trisca,t,1)
202  do k=k1,k2
203  dati%i = levs + 2 + (k-1) * 2 + 1 ! Divergence
204  dati%f => trisca(:,k)
205  call sigio_rrdati(lusig,head,dati,irets)
206  if(irets /= 0) return
207  dati%i = levs + 2 + (k-1) * 2 + 2 ! Vorticity
208  dati%f => triscb(:,k)
209  call sigio_rrdati(lusig,head,dati,irets)
210  if(irets /= 0) return
211  enddo
212  call sptezmv(0,jcap,idrt,io,jo,klen,trisca(1,k1),triscb(1,k1), &
213  u(1,k1),v(1,k1),1)
214  call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),d(1,k1),1)
215 
216 ! call sptezm(0,jcap,idrt,io,jo,1,triscb,z(1,k),1)
217  write(*,*)' retriving d/z for k=',k,' k1=',k1,' k2=',k2
218 ! datm%z(3,:) = datm%z(3,:)+2*con_omega/sqrt(1.5)
219 ! call sptezm(0,jcap,idrt,io,jo,klen,datm%z,z,1)
220  write(*,*)' start get tracer'
221  do nt=1,ntrac
222  do k=k1,k2
223  dati%i = levs * (2+nt) + 2 + k ! Tracers starting with q
224  dati%f => trisca(:,k)
225  call sigio_rrdati(lusig,head,dati,irets)
226  enddo
227  call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),trc(1,k1,nt),1)
228  write(*,*)' retriving d/z for nt=',nt,'ntrac=',ntrac,'k=',k,' k1=',k1,' k2=',k2
229  enddo
230  !t=t/(1+con_fvirt*sh)
231  write(*,*)' end get tracer,idvm=',idvm,'ijo=',ijo,'ntrac=',ntrac
232 !
233 !-- get temp
234  if (mod(idvm/10,10) == 3) then ! Enthalpy case
235  allocate(cpi(0:ntrac))
236 ! write(*,*)'aft read sig, cpi=',head%cpi
237  cpi(0:ntrac) = head%cpi(1:ntrac+1)
238 ! write(*,*)'cpi=',cpi(0:ntrac)
239 !$omp parallel do private(k,i,xcp,sumq,n)
240  do k=k1,k2
241  do i=1,ijo
242  xcp = 0.0
243  sumq = 0.0
244  do n=1,ntrac
245  if( cpi(n) /= 0.0 ) then
246  xcp = xcp + cpi(n)*trc(i,k,n)
247  sumq = sumq + trc(i,k,n)
248  endif
249  enddo
250  xcp = (1.-sumq)*cpi(0) + xcp
251  t(i,k) = t(i,k) / xcp ! Now g1 contains T
252  enddo
253  enddo
254  if (allocated(cpi)) deallocate(cpi)
255  else
256 !$omp parallel do private(i,k)
257  do k=k1,k2
258  do i=1,ijo
259  t(i,k) = t(i,k) / (1+con_fvirt*trc(i,k,1)) !get temp from virtual temp
260  enddo
261  enddo
262  endif
263  endif
264 ! write(*,*)'end comput t'
265  iret=0
266 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267 end subroutine
268 
269 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
297  subroutine modstuff(km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,&
298  pi,pm,om)
299  use sigio_module, only: sigio_modprd
300  implicit none
301  integer,intent(in):: km,idvc,idsl,nvcoord
302  real,intent(in):: vcoord(km+1,nvcoord)
303  real,intent(in):: ps,psx,psy
304  real,intent(in):: u(km),v(km),d(km)
305 ! real,intent(out):: pi(km+1),pm(km)
306  real*8, intent(out):: pi(km+1),pm(km)
307  real,intent(out):: om(km)
308  real*8 ps8,pm8(km),pd8(km),vcoord8(km+1,nvcoord)
309  real*8 dpmdps(km),dpddps(km),dpidps(km+1),pi8(km+1)
310  real vgradp,pd(km),px(km),py(km),os
311  integer k,iret,logk
312 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
313  ps8=ps
314  vcoord8=vcoord
315  call sigio_modprd(1,1,km,nvcoord,idvc,idsl,vcoord8,iret,&
316  ps=(/ps8/),&
317  pm=pm8,pd=pd8,dpmdps=dpmdps,dpddps=dpddps)
318 !
319 !jw: has to be 8 real for wam
320  pi8(1)=ps
321  pm=pm8
322 ! pd=pd8
323  dpidps(1)=1.
324  do k=1,km
325  pi8(k+1)=pi8(k)-pd8(k)
326  dpidps(k+1)=dpidps(k)-dpddps(k)
327 ! if(pi(8)<0.) then
328 ! print *,'in modstuff,pi8=',pi8(k)
329 ! endif
330  enddo
331  pi=pi8
332 !
333  os=0
334  do k=km,1,-1
335  vgradp=u(k)*psx+v(k)*psy
336  os=os-vgradp*ps*(dpmdps(k)-dpidps(k+1))-d(k)*(pm(k)-pi(k+1))
337  om(k)=vgradp*ps*dpmdps(k)+os
338  os=os-vgradp*ps*(dpidps(k)-dpmdps(k))-d(k)*(pi(k)-pm(k))
339  enddo
340  px=ps*dpmdps*psx
341  py=ps*dpmdps*psy
342  end subroutine
343 
344 !-------------------------------------------------------------------------------
376  subroutine modstuff2(im,ix,km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,&
377  pi,pm,om,me)
378  use sigio_module, only : sigio_modprd
379  implicit none
380  integer, intent(in) :: im,ix,km,idvc,idsl,nvcoord,me
381  real, intent(in) :: vcoord(km+1,nvcoord)
382  real, dimension(ix), intent(in) :: ps,psx,psy
383  real, dimension(ix,km), intent(in) :: u,v,d
384  real*8, dimension(ix,km+1), intent(out) :: pi
385  real*8, dimension(ix,km), intent(out) :: pm
386  real, dimension(ix,km), intent(out) :: om
387 ! real*8, allocatable :: ps8(:), pm8(:,:), pd8(:,:),dpmdps(:,:),dpddps(:,:), &
388 ! dpidps(:,:),pi8(:,:),vcoord8(:,:)
389 ! real, allocatable :: os(:)
390 ! real, allocatable :: pd(:,:),px(:,:), py(:,:), os(:)
391 
392 ! real vgradpps
393 
394  real*8 ps8(ix),pm8(ix,km),pd8(ix,km),vcoord8(km+1,nvcoord)
395  real*8 dpmdps(ix,km),dpddps(ix,km),dpidps(ix,km+1),pi8(ix,km+1)
396  real vgradpps,pd(im,km),os(im)
397 ! real vgradpps,pd(im,km),px(im,km),py(im,km),os(im),tem
398  integer i,k,iret,logk
399 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
400 
401  ps8 = ps
402  vcoord8 = vcoord
403  call sigio_modprd(im,ix,km,nvcoord,idvc,idsl,vcoord8,iret, &
404  ps=ps8,pd=pd8,dpddps=dpddps,pm=pm8,dpmdps=dpmdps)
405 
406 !
407 ! if (me == 0) then
408 ! write(*,*)' pd8=',pd8(1,60:64)
409 ! write(*,*)' pm8=',pm8(1,60:64)
410 ! endif
411 !jw: has to be 8 real for wam
412 
413 !$omp parallel do private(i)
414  do i=1,im
415  pi8(i,1) = ps(i)
416  dpidps(i,1) = 1.
417  os(i) = 0
418  pi(i,1) = pi8(i,1)
419  enddo
420  do k=1,km
421 !$omp parallel do private(i)
422  do i=1,im
423  pi8(i,k+1) = pi8(i,k) - pd8(i,k)
424  dpidps(i,k+1) = dpidps(i,k) - dpddps(i,k)
425 ! if(pi(i,8)<0.) then
426 ! print *,'in modstuff,pi8=',pi8(i,k),' i=',i,' k=',k,' me=',me
427 ! endif
428  pi(i,k+1) = pi8(i,k+1)
429  pm(i,k) = pm8(i,k)
430  enddo
431  enddo
432 !
433  do k=km,1,-1
434 !$omp parallel do private(i,vgradpps)
435  do i=1,im
436  vgradpps = (u(i,k)*psx(i) + v(i,k)*psy(i)) * ps(i)
437 
438  os(i) = os(i) - vgradpps*(dpmdps(i,k)-dpidps(i,k+1)) &
439  - d(i,k)*(pm(i,k)-pi(i,k+1))
440 
441  om(i,k) = os(i) + vgradpps*dpmdps(i,k)
442 
443  os(i) = os(i) - vgradpps*(dpidps(i,k)-dpmdps(i,k)) &
444  - d(i,k)*(pi(i,k)-pm(i,k))
445  enddo
446  enddo
447  end subroutine
448 
449 !-----------------------------------------------------------------------
497  subroutine trssc(jcap,nc,km,ntrac,idvc,idvm,idsl,nvcoord,vcoord, &
498  cpi,idrt,lonb,latb,ijl,ijn,j1,j2,jc,chgq0, &
499  szs,sps,st,sd,sz,sq,gfszs,gfsps,gfsp,gfsdp, &
500  gfst,gfsu,gfsv,gfsq,gfsw)
501  implicit none
502  integer,intent(in)::jcap,nc,km,ntrac,idvc,idvm,idsl,nvcoord,idrt,lonb,latb
503  integer,intent(in)::ijl,ijn,j1,j2,jc,chgq0
504  real,intent(in):: szs(nc),sps(nc),st(nc,km),sd(nc,km),sz(nc,km),sq(nc,km*ntrac)
505  real,intent(in):: cpi(0:ntrac)
506  real*8,intent(in):: vcoord(km+1,nvcoord)
507  real,dimension(ijn),intent(inout):: gfszs,gfsps
508  real,dimension(ijn,km),intent(inout):: gfsp,gfsdp,gfst,gfsu,gfsv,gfsw
509  real,dimension(ijn,km*ntrac),intent(inout):: gfsq
510  real zs(ijl),ps(ijl),t(ijl,km),u(ijl,km),v(ijl,km),q(ijl,km*ntrac)
511  real wi(ijl,km),pi(ijl,km),dpo(ijl,km)
512  real tvcon,xcp,sumq
513  integer thermodyn_id,jn,js,is,in
514  integer jj,jjm,k,n,j,i,ij,lonb2
515 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
516 ! spectral transforms
517  if(j1==732)print*,'sample input to trssc= ',jcap,nc,km,ntrac, &
518  idvc,idvm,idsl,nvcoord, &
519  idrt,lonb,latb,ijl,ijn,j1,j2,jc,chgq0
520  lonb2=lonb*2
521  ij=lonb2*(j2-j1+1)
522  in=1
523  is=1+lonb
524  call sptran(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijl, &
525  j1,j2,jc,szs,zs(in),zs(is),1)
526  call sptran(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijl, &
527  j1,j2,jc,sps,ps(in),ps(is),1)
528  call sptran(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijl, &
529  j1,j2,jc,st,t(in,1),t(is,1),1)
530  call sptranv(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijl, &
531  j1,j2,jc,sd,sz,u(in,1),u(is,1),v(in,1),v(is,1),1)
532  call sptran(0,jcap,idrt,lonb,latb,km*ntrac,1,1,lonb2,lonb2,nc,ijl, &
533  j1,j2,jc,sq,q(in,1),q(is,1),1)
534  if(j1==732)then
535  do k=1,km
536  do i=1,ijl
537  if(t(i,k)>400. .or. t(i,k)<100.)print*,'bad T from sptran',i,k,t(i,k)
538  if(q(i,k)>1.)print*,'bad Q from sptran',i,k,q(i,k)
539  if(q(i,2*k)>1.)print*,'bad Q from sptran',i,k,q(i,2*k)
540  if(q(i,3*k)>1.)print*,'bad Q from sptran',i,k,q(i,3*k)
541  end do
542  end do
543  end if
544  select case(mod(idvm,10))
545  case(0,1)
546  do i=1,ij
547  ps(i)=1.e3*exp(ps(i))
548  enddo
549  case(2)
550  do i=1,ij
551  ps(i)=1.e3*ps(i)
552  enddo
553  case default
554  do i=1,ij
555  ps(i)=1.e3*exp(ps(i))
556  enddo
557  end select
558 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
559  thermodyn_id=mod(idvm/10,10)
560  if (thermodyn_id == 3) then
561  do k=1,km
562  do i=1,ij
563  t(i,k) = t(i,k)/cpi(0) ! enthalpy (cpt/cpd)
564  end do
565  end do
566 !
567  endif
568 
569  call getomega(jcap,nc,km,idvc,idvm,idrt,idsl, &
570  nvcoord,vcoord,lonb,latb,ijl,j1,j2,1,sd,sps, &
571  ps,t,u,v,wi,pi,dpo)
572  if(j1==732)then
573  do k=1,km
574  do i=1,ijl
575  if(t(i,k)>400. .or. t(i,k)<100.)print*,'bad T after getomega',i,k,t(i,k)
576  if(q(i,k)>1. )print*,'bad Q after getomega',i,k,q(i,k)
577  if(q(i,2*k)>1. )print*,'bad Q after getomega',i,2*k,q(i,2*k)
578  end do
579  end do
580  end if
581  if(thermodyn_id /= 2)then
582 ! convert to surface pressure and temperature
583  if (thermodyn_id == 3) then
584  do k=1,km
585  do i=1,ij
586  xcp = 0.0
587  sumq = 0.0
588  do n=1,ntrac
589  if( cpi(n) /= 0.0 ) then
590  xcp = xcp + cpi(n)*q(i,k+(n-1)*km)
591  sumq = sumq + q(i,k+(n-1)*km)
592  endif
593  enddo
594  t(i,k) = t(i,k)/((1.-sumq)*cpi(0)+xcp)
595  end do
596  end do
597 
598  else
599  tvcon=(461.50/287.05-1.)
600  t(:,:) = t(:,:)/(1.+tvcon*q(:,1:km))
601  endif
602  end if
603  if(j1==732)then
604  do k=1,km
605  do i=1,ijl
606  if(t(i,k)>400. .or. t(i,k)<100.)print*,'bad T after Tv to T',i,k,t(i,k)
607  if(q(i,k)>1.)print*,'bad Q after Tv to T',i,k,q(i,k)
608  if(q(i,2*k)>1. )print*,'bad Q after Tv to T',i,k,q(i,2*k)
609  end do
610  end do
611  end if
612 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
613 !----force tracers to be positive
614  if (chgq0 == 1) q = max(q, 0.0)
615 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
616 ! pass data to gfsdatao
617  do j=1,j2-j1+1
618  jn=j+j1-1
619  js=latb+1-jn
620  jn=(jn-1)*lonb
621  js=(js-1)*lonb
622  jj=j*lonb
623  jjm=(j-1)*lonb
624  do i=1,lonb
625  gfszs(i+jn) = zs(i+jjm)
626  gfsps(i+jn) = ps(i+jjm)
627  enddo
628  do i=1,lonb
629  gfszs(i+js) = zs(i+jj)
630  gfsps(i+js) = ps(i+jj)
631  enddo
632  do k=1,km
633  do i=1,lonb
634  gfsdp(i+jn,k) = dpo(i+jjm,k)
635  gfsp(i+jn,k) = pi(i+jjm,k)
636  gfst(i+jn,k) = t(i+jjm,k)
637  gfsu(i+jn,k) = u(i+jjm,k)
638  gfsv(i+jn,k) = v(i+jjm,k)
639  gfsw(i+jn,k) = wi(i+jjm,k)
640  enddo
641  do i=1,lonb
642  gfsdp(i+js,k) = dpo(i+jj,k)
643  gfsp(i+js,k) = pi(i+jj,k)
644  gfst(i+js,k) = t(i+jj,k)
645  gfsu(i+js,k) = u(i+jj,k)
646  gfsv(i+js,k) = v(i+jj,k)
647  gfsw(i+js,k) = wi(i+jj,k)
648  enddo
649  enddo
650  do k=1,km*ntrac
651  do i=1,lonb
652  gfsq(i+jn,k) = q(i+jjm,k)
653  enddo
654  do i=1,lonb
655  gfsq(i+js,k) = q(i+jj,k)
656  enddo
657  enddo
658  enddo
659  return
660  end
661 !-----------------------------------------------------------------------
689  subroutine getomega(jcap,nc,km,idvc,idvm,idrt,idsl,nvcoord,vcoord, &
690  lonb,latb,ijn,j1,j2,jc,sd,sps,psi,ti,ui,vi,wi,pm,pd)
691 !!!!!
692  use sigio_module, only : sigio_modprd
693  implicit none
694 
695  integer,intent(in):: jcap,nc,km,idvc,idvm,idrt,idsl,nvcoord
696  integer,intent(in):: lonb,latb,j1,j2,jc,ijn
697  real*8,intent(in):: vcoord(km+1,nvcoord)
698  real,intent(in):: sd(nc,km),sps(nc)
699  real,intent(in):: psi(ijn),ti(ijn,km),ui(ijn,km),vi(ijn,km)
700  real,intent(out):: wi(ijn,km),pm(ijn,km),pd(ijn,km)
701  real :: pi(ijn,km+1)
702  real :: os
703  real*8 psi8(ijn),ti8(ijn,km),pm8(ijn,km),pd8(ijn,km)
704  real*8 dpmdps(ijn,km),dpddps(ijn,km),dpidps(ijn,km+1),vgradp,psmean
705  real di(ijn,km),psx(ijn),psy(ijn)
706  integer k,i,ij,lonb2,iret,is,in
707 !----1. spectral transform
708  lonb2=lonb*2
709  ij=lonb2*(j2-j1+1)
710  in=1
711  is=1+lonb
712  call sptrand(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijn, &
713  j1,j2,jc,sps,psmean,psx(in),psx(is),psy(in),psy(is),1)
714 
715  call sptran(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijn, &
716  j1,j2,jc,sd,di(in,1),di(is,1),1)
717  psi8=psi
718  ti8=ti
719 
720  call sigio_modprd(ijn,ijn,km,nvcoord,idvc,idsl,vcoord,iret, &
721  ps=psi8,t=ti8,pm=pm8,pd=pd8,dpmdps=dpmdps,dpddps=dpddps)
722  pm=pm8
723  pd=pd8
724 
725  select case(mod(idvm,10))
726  case(0,1)
727  continue
728  case(2)
729  do i=1,ijn
730  psx(i)=psx(i)/(psi(i)*1.0e-3)
731  psy(i)=psy(i)/(psi(i)*1.0e-3)
732  enddo
733  case default
734  do i=1,ijn
735  psx(i)=psx(i)/psi(i)
736  psy(i)=psy(i)/psi(i)
737  enddo
738  end select
739 
740 !----3.omeda from modstuff
741  do i=1,ijn
742  pi(i,1)=psi(i)
743  dpidps(i,1)=1.
744  enddo
745  do k=1,km
746  do i=1,ijn
747  pi(i,k+1)=pi(i,k)-pd(i,k)
748  dpidps(i,k+1)=dpidps(i,k)-dpddps(i,k)
749  enddo
750  enddo
751  do i=1,ijn
752  os=0.
753  do k=km,1,-1
754  vgradp=ui(i,k)*psx(i)+vi(i,k)*psy(i)
755  os=os-vgradp*psi(i)*(dpmdps(i,k)-dpidps(i,k+1))- &
756  di(i,k)*(pm(i,k)-pi(i,k+1))
757  wi(i,k)=vgradp*psi(i)*dpmdps(i,k)+os
758  os=os-vgradp*psi(i)*(dpidps(i,k)-dpmdps(i,k))- &
759  di(i,k)*(pi(i,k)-pm(i,k))
760  enddo
761 !
762  enddo
763 !---
764  return
765  end subroutine
rtsig
subroutine rtsig(lusig, head, k1, k2, kgds, ijo, levs, ntrac, jcap, lnt2, me, h, p, px, py, t, u, v, d, trc, iret)
rtsig() reads and transforms sigma file.
Definition: GFSPOSTSIG.F:94
physcons_post
Definition: physcons.f:1
modstuff2
subroutine modstuff2(im, ix, km, idvc, idsl, nvcoord, vcoord, ps, psx, psy, d, u, v, pi, pm, om, me)
modstuff2() computes model coordinate dependent functions.
Definition: GFSPOSTSIG.F:376
getomega
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)
getomega() computes omega.
Definition: GFSPOSTSIG.F:689
modstuff
subroutine modstuff(km, idvc, idsl, nvcoord, vcoord, ps, psx, psy, d, u, v, pi, pm, om)
modstuff() computes model coordinate dependent functions.
Definition: GFSPOSTSIG.F:297
trssc
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)
trssc() transforms sigma spectral fields to grid.
Definition: GFSPOSTSIG.F:497