UPP  001
 All Data Structures Files Functions Pages
GFSPOSTSIG.F
Go to the documentation of this file.
1 
52 subroutine rtsig(lusig,head,k1,k2,kgds,ijo,levs,ntrac,jcap,lnt2,me, &
53  h,p,px,py,t,u,v,d,trc,iret)
54 
55  use sigio_module, only : sigio_intkind, sigio_head
56  use sigio_r_module, only : sigio_dati, sigio_rrdati
57  use physcons_post, only : con_omega, con_fvirt
58  use omp_lib
59  implicit none
60  integer(sigio_intkind),intent(in) :: lusig
61  type(sigio_head), intent(in) :: head
62  integer,intent(in) :: k1,k2,kgds(200),ijo,levs,ntrac,jcap,lnt2,me
63  real,dimension(ijo), intent(out) :: h,p,px,py
64  real,dimension(ijo,k1:k2),intent(out):: t,u,v,d
65  real,dimension(ijo,k1:k2,ntrac),intent(out),target :: trc
66  integer,intent(out) :: iret
67 !
68  integer idrt,io,jo,iidea
69 ! integer idrt,io,jo,mo3,mct,iidea,midea
70  integer(sigio_intkind):: irets
71 ! type(sigio_datm):: datm
72  type(sigio_dati) dati
73 ! type griddata
74 ! real,dimension(:,:),pointer :: datm
75 ! endtype griddata
76 ! type(griddata),dimension(:),pointer :: datatrc
77  real, target :: trisca(lnt2,k1:k2+1), triscb(lnt2,k1:k2)
78  real,dimension(:), allocatable :: cpi
79  real,dimension(:,:),allocatable :: wrk
80  integer io3,ict,jct,n,i,k,jc,nt
81  integer idvm, klen
82  real pmean,sumq,xcp
83 ! integer, parameter :: latch=20
84 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85 ! Determine output grid
86  idrt = kgds(1)
87  if(kgds(1) == 0 .and. kgds(4) < 90000) idrt = 256
88  io = kgds(2)
89  jo = kgds(3)
90 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91 ! Read and transform surface fields
92  iret = 1
93  if (me == 0) then
94  print*,'Debug rtsig= ',lusig,k1,k2,ijo,kgds(1:20)
95  endif
96 
97  idvm = head%idvm
98 ! jc = omp_get_num_threads()
99 ! write(0,*)' in RTSIG lnt2=',lnt2,' threads=',jc,' latch=',latch, &
100 ! ' jcap=',jcap,' io=',io,' jo=',jo,' ijo=',ijo
101 !
102  if (k2 < k1) return
103 
104  dati%i = 1 ! hs
105  dati%f => trisca(:,k1)
106  call sigio_rrdati(lusig,head,dati,irets)
107  if(irets /= 0) return
108 
109 ! call sptez(0,jcap,idrt,io,jo,trisca(1,k1),h,1)
110 ! call sptez(0,jcap,idrt,io,jo,dats%hs,h,1)
111 ! call sptez(0,jcap,idrt,io,jo,dats%ps,p,1)
112 ! call sptezj(jcap,lnt2,1,idrt,io,jo,jc,trisca,h,latch,1)
113 !
114  dati%i = 2 ! Surface pressure
115  dati%f => trisca(:,k1+1)
116  call sigio_rrdati(lusig,head,dati,irets)
117  if(irets /= 0) return
118 !
119 ! call sptez(0,jcap,idrt,io,jo,trisca(1,k1),p,1)
120 ! call sptezj(jcap,lnt2,1,idrt,io,jo,jc,trisca,p,latch,1)
121 !--
122  allocate(wrk(ijo,2))
123  call sptezm(0,jcap,idrt,io,jo,2,trisca(1,k1),wrk,1)
124  if( mod(idvm,10) < 2) then
125 !$omp parallel do private(i)
126  do i=1,ijo
127  h(i) = wrk(i,1)
128  p(i) = 1.e3*exp(wrk(i,2))
129 ! p(i) = 1.e3*exp(p(i))
130  enddo
131  elseif(mod(idvm,10) == 2) then
132 !$omp parallel do private(i)
133  do i=1,ijo
134  h(i) = wrk(i,1)
135  p(i) = 1000.*wrk(i,2)
136 ! p(i) = 1000.*p(i)
137  enddo
138  endif
139  if (allocated(wrk)) deallocate(wrk)
140 
141  call sptezd(0,jcap,idrt,io,jo,trisca(1,k1+1),pmean,px,py,1)
142  iret = 0
143 
144 ! if (k2 < k1) return
145 
146 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147 ! Read and transform fields on levels k1 through k2
148  iret = 2
149  if (k2 >= k1) then
150  klen = k2-k1+1
151  do k=k1,k2
152  write(0,*)' retriving T for k=',k,' k1=',k1,' k2=',k2
153  dati%i = k + 2 ! Virtual Temperature or CpT
154  dati%f => trisca(:,k)
155 
156  call sigio_rrdati(lusig,head,dati,iret)
157  enddo
158  call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),t(1,k1),1)
159 ! call sptezm(0,jcap,idrt,io,jo,klen,trisca,t,1)
160  do k=k1,k2
161  dati%i = levs + 2 + (k-1) * 2 + 1 ! Divergence
162  dati%f => trisca(:,k)
163  call sigio_rrdati(lusig,head,dati,irets)
164  if(irets /= 0) return
165  dati%i = levs + 2 + (k-1) * 2 + 2 ! Vorticity
166  dati%f => triscb(:,k)
167  call sigio_rrdati(lusig,head,dati,irets)
168  if(irets /= 0) return
169  enddo
170  call sptezmv(0,jcap,idrt,io,jo,klen,trisca(1,k1),triscb(1,k1), &
171  u(1,k1),v(1,k1),1)
172  call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),d(1,k1),1)
173 
174 ! call sptezm(0,jcap,idrt,io,jo,1,triscb,z(1,k),1)
175  write(0,*)' retriving d/z for k=',k,' k1=',k1,' k2=',k2
176 ! datm%z(3,:) = datm%z(3,:)+2*con_omega/sqrt(1.5)
177 ! call sptezm(0,jcap,idrt,io,jo,klen,datm%z,z,1)
178  write(0,*)' start get tracer'
179  do nt=1,ntrac
180  do k=k1,k2
181  dati%i = levs * (2+nt) + 2 + k ! Tracers starting with q
182  dati%f => trisca(:,k)
183  call sigio_rrdati(lusig,head,dati,irets)
184  enddo
185  call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),trc(1,k1,nt),1)
186  write(0,*)' retriving d/z for nt=',nt,'ntrac=',ntrac,'k=',k,' k1=',k1,' k2=',k2
187  enddo
188  !t=t/(1+con_fvirt*sh)
189  write(0,*)' end get tracer,idvm=',idvm,'ijo=',ijo,'ntrac=',ntrac
190 !
191 !-- get temp
192  if (mod(idvm/10,10) == 3) then ! Enthalpy case
193  allocate(cpi(0:ntrac))
194 ! write(0,*)'aft read sig, cpi=',head%cpi
195  cpi(0:ntrac) = head%cpi(1:ntrac+1)
196 ! write(0,*)'cpi=',cpi(0:ntrac)
197 !$omp parallel do private(k,i,xcp,sumq,n)
198  do k=k1,k2
199  do i=1,ijo
200  xcp = 0.0
201  sumq = 0.0
202  do n=1,ntrac
203  if( cpi(n) /= 0.0 ) then
204  xcp = xcp + cpi(n)*trc(i,k,n)
205  sumq = sumq + trc(i,k,n)
206  endif
207  enddo
208  xcp = (1.-sumq)*cpi(0) + xcp
209  t(i,k) = t(i,k) / xcp ! Now g1 contains T
210  enddo
211  enddo
212  if (allocated(cpi)) deallocate(cpi)
213  else
214 !$omp parallel do private(i,k)
215  do k=k1,k2
216  do i=1,ijo
217  t(i,k) = t(i,k) / (1+con_fvirt*trc(i,k,1)) !get temp from virtual temp
218  enddo
219  enddo
220  endif
221  endif
222 ! write(0,*)'end comput t'
223  iret=0
224 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225 end subroutine
226 
227 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255  subroutine modstuff(km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,&
256  pi,pm,om)
257  use sigio_module, only: sigio_modprd
258  implicit none
259  integer,intent(in):: km,idvc,idsl,nvcoord
260  real,intent(in):: vcoord(km+1,nvcoord)
261  real,intent(in):: ps,psx,psy
262  real,intent(in):: u(km),v(km),d(km)
263 ! real,intent(out):: pi(km+1),pm(km)
264  real*8, intent(out):: pi(km+1),pm(km)
265  real,intent(out):: om(km)
266  real*8 ps8,pm8(km),pd8(km),vcoord8(km+1,nvcoord)
267  real*8 dpmdps(km),dpddps(km),dpidps(km+1),pi8(km+1)
268  real vgradp,pd(km),px(km),py(km),os
269  integer k,iret,logk
270 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
271  ps8=ps
272  vcoord8=vcoord
273  call sigio_modprd(1,1,km,nvcoord,idvc,idsl,vcoord8,iret,&
274  ps=(/ps8/),&
275  pm=pm8,pd=pd8,dpmdps=dpmdps,dpddps=dpddps)
276 !
277 !jw: has to be 8 real for wam
278  pi8(1)=ps
279  pm=pm8
280 ! pd=pd8
281  dpidps(1)=1.
282  do k=1,km
283  pi8(k+1)=pi8(k)-pd8(k)
284  dpidps(k+1)=dpidps(k)-dpddps(k)
285 ! if(pi(8)<0.) then
286 ! print *,'in modstuff,pi8=',pi8(k)
287 ! endif
288  enddo
289  pi=pi8
290 !
291  os=0
292  do k=km,1,-1
293  vgradp=u(k)*psx+v(k)*psy
294  os=os-vgradp*ps*(dpmdps(k)-dpidps(k+1))-d(k)*(pm(k)-pi(k+1))
295  om(k)=vgradp*ps*dpmdps(k)+os
296  os=os-vgradp*ps*(dpidps(k)-dpmdps(k))-d(k)*(pi(k)-pm(k))
297  enddo
298  px=ps*dpmdps*psx
299  py=ps*dpmdps*psy
300  end subroutine
301 
302 !-------------------------------------------------------------------------------
333  subroutine modstuff2(im,ix,km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,&
334  pi,pm,om,me)
335  use sigio_module, only : sigio_modprd
336  implicit none
337  integer, intent(in) :: im,ix,km,idvc,idsl,nvcoord,me
338  real, intent(in) :: vcoord(km+1,nvcoord)
339  real, dimension(ix), intent(in) :: ps,psx,psy
340  real, dimension(ix,km), intent(in) :: u,v,d
341  real*8, dimension(ix,km+1), intent(out) :: pi
342  real*8, dimension(ix,km), intent(out) :: pm
343  real, dimension(ix,km), intent(out) :: om
344 ! real*8, allocatable :: ps8(:), pm8(:,:), pd8(:,:),dpmdps(:,:),dpddps(:,:), &
345 ! dpidps(:,:),pi8(:,:),vcoord8(:,:)
346 ! real, allocatable :: os(:)
347 ! real, allocatable :: pd(:,:),px(:,:), py(:,:), os(:)
348 
349 ! real vgradpps
350 
351  real*8 ps8(ix),pm8(ix,km),pd8(ix,km),vcoord8(km+1,nvcoord)
352  real*8 dpmdps(ix,km),dpddps(ix,km),dpidps(ix,km+1),pi8(ix,km+1)
353  real vgradpps,pd(im,km),os(im)
354 ! real vgradpps,pd(im,km),px(im,km),py(im,km),os(im),tem
355  integer i,k,iret,logk
356 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
357 
358  ps8 = ps
359  vcoord8 = vcoord
360  call sigio_modprd(im,ix,km,nvcoord,idvc,idsl,vcoord8,iret, &
361  ps=ps8,pd=pd8,dpddps=dpddps,pm=pm8,dpmdps=dpmdps)
362 
363 !
364 ! if (me == 0) then
365 ! write(0,*)' pd8=',pd8(1,60:64)
366 ! write(0,*)' pm8=',pm8(1,60:64)
367 ! endif
368 !jw: has to be 8 real for wam
369 
370 !$omp parallel do private(i)
371  do i=1,im
372  pi8(i,1) = ps(i)
373  dpidps(i,1) = 1.
374  os(i) = 0
375  pi(i,1) = pi8(i,1)
376  enddo
377  do k=1,km
378 !$omp parallel do private(i)
379  do i=1,im
380  pi8(i,k+1) = pi8(i,k) - pd8(i,k)
381  dpidps(i,k+1) = dpidps(i,k) - dpddps(i,k)
382 ! if(pi(i,8)<0.) then
383 ! print *,'in modstuff,pi8=',pi8(i,k),' i=',i,' k=',k,' me=',me
384 ! endif
385  pi(i,k+1) = pi8(i,k+1)
386  pm(i,k) = pm8(i,k)
387  enddo
388  enddo
389 !
390  do k=km,1,-1
391 !$omp parallel do private(i,vgradpps)
392  do i=1,im
393  vgradpps = (u(i,k)*psx(i) + v(i,k)*psy(i)) * ps(i)
394 
395  os(i) = os(i) - vgradpps*(dpmdps(i,k)-dpidps(i,k+1)) &
396  - d(i,k)*(pm(i,k)-pi(i,k+1))
397 
398  om(i,k) = os(i) + vgradpps*dpmdps(i,k)
399 
400  os(i) = os(i) - vgradpps*(dpidps(i,k)-dpmdps(i,k)) &
401  - d(i,k)*(pi(i,k)-pm(i,k))
402  enddo
403  enddo
404  end subroutine
405 
406 !-----------------------------------------------------------------------
444  subroutine trssc(jcap,nc,km,ntrac,idvc,idvm,idsl,nvcoord,vcoord, &
445  cpi,idrt,lonb,latb,ijl,ijn,j1,j2,jc,chgq0, &
446  szs,sps,st,sd,sz,sq,gfszs,gfsps,gfsp,gfsdp, &
447  gfst,gfsu,gfsv,gfsq,gfsw)
448  implicit none
449  integer,intent(in)::jcap,nc,km,ntrac,idvc,idvm,idsl,nvcoord,idrt,lonb,latb
450  integer,intent(in)::ijl,ijn,j1,j2,jc,chgq0
451  real,intent(in):: szs(nc),sps(nc),st(nc,km),sd(nc,km),sz(nc,km),sq(nc,km*ntrac)
452  real,intent(in):: cpi(0:ntrac)
453  real*8,intent(in):: vcoord(km+1,nvcoord)
454  real,dimension(ijn),intent(inout):: gfszs,gfsps
455  real,dimension(ijn,km),intent(inout):: gfsp,gfsdp,gfst,gfsu,gfsv,gfsw
456  real,dimension(ijn,km*ntrac),intent(inout):: gfsq
457  real zs(ijl),ps(ijl),t(ijl,km),u(ijl,km),v(ijl,km),q(ijl,km*ntrac)
458  real wi(ijl,km),pi(ijl,km),dpo(ijl,km)
459  real tvcon,xcp,sumq
460  integer thermodyn_id,jn,js,is,in
461  integer jj,jjm,k,n,j,i,ij,lonb2
462 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
463 ! spectral transforms
464  if(j1==732)print*,'sample input to trssc= ',jcap,nc,km,ntrac, &
465  idvc,idvm,idsl,nvcoord, &
466  idrt,lonb,latb,ijl,ijn,j1,j2,jc,chgq0
467  lonb2=lonb*2
468  ij=lonb2*(j2-j1+1)
469  in=1
470  is=1+lonb
471  call sptran(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijl, &
472  j1,j2,jc,szs,zs(in),zs(is),1)
473  call sptran(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijl, &
474  j1,j2,jc,sps,ps(in),ps(is),1)
475  call sptran(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijl, &
476  j1,j2,jc,st,t(in,1),t(is,1),1)
477  call sptranv(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijl, &
478  j1,j2,jc,sd,sz,u(in,1),u(is,1),v(in,1),v(is,1),1)
479  call sptran(0,jcap,idrt,lonb,latb,km*ntrac,1,1,lonb2,lonb2,nc,ijl, &
480  j1,j2,jc,sq,q(in,1),q(is,1),1)
481  if(j1==732)then
482  do k=1,km
483  do i=1,ijl
484  if(t(i,k)>400. .or. t(i,k)<100.)print*,'bad T from sptran',i,k,t(i,k)
485  if(q(i,k)>1.)print*,'bad Q from sptran',i,k,q(i,k)
486  if(q(i,2*k)>1.)print*,'bad Q from sptran',i,k,q(i,2*k)
487  if(q(i,3*k)>1.)print*,'bad Q from sptran',i,k,q(i,3*k)
488  end do
489  end do
490  end if
491  select case(mod(idvm,10))
492  case(0,1)
493  do i=1,ij
494  ps(i)=1.e3*exp(ps(i))
495  enddo
496  case(2)
497  do i=1,ij
498  ps(i)=1.e3*ps(i)
499  enddo
500  case default
501  do i=1,ij
502  ps(i)=1.e3*exp(ps(i))
503  enddo
504  end select
505 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
506  thermodyn_id=mod(idvm/10,10)
507  if (thermodyn_id == 3) then
508  do k=1,km
509  do i=1,ij
510  t(i,k) = t(i,k)/cpi(0) ! enthalpy (cpt/cpd)
511  end do
512  end do
513 !
514  endif
515 
516  call getomega(jcap,nc,km,idvc,idvm,idrt,idsl, &
517  nvcoord,vcoord,lonb,latb,ijl,j1,j2,1,sd,sps, &
518  ps,t,u,v,wi,pi,dpo)
519  if(j1==732)then
520  do k=1,km
521  do i=1,ijl
522  if(t(i,k)>400. .or. t(i,k)<100.)print*,'bad T after getomega',i,k,t(i,k)
523  if(q(i,k)>1. )print*,'bad Q after getomega',i,k,q(i,k)
524  if(q(i,2*k)>1. )print*,'bad Q after getomega',i,2*k,q(i,2*k)
525  end do
526  end do
527  end if
528  if(thermodyn_id /= 2)then
529 ! convert to surface pressure and temperature
530  if (thermodyn_id == 3) then
531  do k=1,km
532  do i=1,ij
533  xcp = 0.0
534  sumq = 0.0
535  do n=1,ntrac
536  if( cpi(n) /= 0.0 ) then
537  xcp = xcp + cpi(n)*q(i,k+(n-1)*km)
538  sumq = sumq + q(i,k+(n-1)*km)
539  endif
540  enddo
541  t(i,k) = t(i,k)/((1.-sumq)*cpi(0)+xcp)
542  end do
543  end do
544 
545  else
546  tvcon=(461.50/287.05-1.)
547  t(:,:) = t(:,:)/(1.+tvcon*q(:,1:km))
548  endif
549  end if
550  if(j1==732)then
551  do k=1,km
552  do i=1,ijl
553  if(t(i,k)>400. .or. t(i,k)<100.)print*,'bad T after Tv to T',i,k,t(i,k)
554  if(q(i,k)>1.)print*,'bad Q after Tv to T',i,k,q(i,k)
555  if(q(i,2*k)>1. )print*,'bad Q after Tv to T',i,k,q(i,2*k)
556  end do
557  end do
558  end if
559 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
560 !----force tracers to be positive
561  if (chgq0 == 1) q = max(q, 0.0)
562 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
563 ! pass data to gfsdatao
564  do j=1,j2-j1+1
565  jn=j+j1-1
566  js=latb+1-jn
567  jn=(jn-1)*lonb
568  js=(js-1)*lonb
569  jj=j*lonb
570  jjm=(j-1)*lonb
571  do i=1,lonb
572  gfszs(i+jn) = zs(i+jjm)
573  gfsps(i+jn) = ps(i+jjm)
574  enddo
575  do i=1,lonb
576  gfszs(i+js) = zs(i+jj)
577  gfsps(i+js) = ps(i+jj)
578  enddo
579  do k=1,km
580  do i=1,lonb
581  gfsdp(i+jn,k) = dpo(i+jjm,k)
582  gfsp(i+jn,k) = pi(i+jjm,k)
583  gfst(i+jn,k) = t(i+jjm,k)
584  gfsu(i+jn,k) = u(i+jjm,k)
585  gfsv(i+jn,k) = v(i+jjm,k)
586  gfsw(i+jn,k) = wi(i+jjm,k)
587  enddo
588  do i=1,lonb
589  gfsdp(i+js,k) = dpo(i+jj,k)
590  gfsp(i+js,k) = pi(i+jj,k)
591  gfst(i+js,k) = t(i+jj,k)
592  gfsu(i+js,k) = u(i+jj,k)
593  gfsv(i+js,k) = v(i+jj,k)
594  gfsw(i+js,k) = wi(i+jj,k)
595  enddo
596  enddo
597  do k=1,km*ntrac
598  do i=1,lonb
599  gfsq(i+jn,k) = q(i+jjm,k)
600  enddo
601  do i=1,lonb
602  gfsq(i+js,k) = q(i+jj,k)
603  enddo
604  enddo
605  enddo
606  return
607  end
608 !-----------------------------------------------------------------------
609  subroutine getomega(jcap,nc,km,idvc,idvm,idrt,idsl,nvcoord,vcoord, &
610  lonb,latb,ijn,j1,j2,jc,sd,sps,psi,ti,ui,vi,wi,pm,pd)
611 !!!!!
612  use sigio_module, only : sigio_modprd
613  implicit none
614 
615  integer,intent(in):: jcap,nc,km,idvc,idvm,idrt,idsl,nvcoord
616  integer,intent(in):: lonb,latb,j1,j2,jc,ijn
617  real*8,intent(in):: vcoord(km+1,nvcoord)
618  real,intent(in):: sd(nc,km),sps(nc)
619  real,intent(in):: psi(ijn),ti(ijn,km),ui(ijn,km),vi(ijn,km)
620  real,intent(out):: wi(ijn,km),pm(ijn,km),pd(ijn,km)
621  real :: pi(ijn,km+1)
622  real :: os
623  real*8 psi8(ijn),ti8(ijn,km),pm8(ijn,km),pd8(ijn,km)
624  real*8 dpmdps(ijn,km),dpddps(ijn,km),dpidps(ijn,km+1),vgradp,psmean
625  real di(ijn,km),psx(ijn),psy(ijn)
626  integer k,i,ij,lonb2,iret,is,in
627 !----1. spectral transform
628  lonb2=lonb*2
629  ij=lonb2*(j2-j1+1)
630  in=1
631  is=1+lonb
632  call sptrand(0,jcap,idrt,lonb,latb,1,1,1,lonb2,lonb2,nc,ijn, &
633  j1,j2,jc,sps,psmean,psx(in),psx(is),psy(in),psy(is),1)
634 
635  call sptran(0,jcap,idrt,lonb,latb,km,1,1,lonb2,lonb2,nc,ijn, &
636  j1,j2,jc,sd,di(in,1),di(is,1),1)
637  psi8=psi
638  ti8=ti
639 
640  call sigio_modprd(ijn,ijn,km,nvcoord,idvc,idsl,vcoord,iret, &
641  ps=psi8,t=ti8,pm=pm8,pd=pd8,dpmdps=dpmdps,dpddps=dpddps)
642  pm=pm8
643  pd=pd8
644 
645  select case(mod(idvm,10))
646  case(0,1)
647  continue
648  case(2)
649  do i=1,ijn
650  psx(i)=psx(i)/(psi(i)*1.0e-3)
651  psy(i)=psy(i)/(psi(i)*1.0e-3)
652  enddo
653  case default
654  do i=1,ijn
655  psx(i)=psx(i)/psi(i)
656  psy(i)=psy(i)/psi(i)
657  enddo
658  end select
659 
660 !----3.omeda from modstuff
661  do i=1,ijn
662  pi(i,1)=psi(i)
663  dpidps(i,1)=1.
664  enddo
665  do k=1,km
666  do i=1,ijn
667  pi(i,k+1)=pi(i,k)-pd(i,k)
668  dpidps(i,k+1)=dpidps(i,k)-dpddps(i,k)
669  enddo
670  enddo
671  do i=1,ijn
672  os=0.
673  do k=km,1,-1
674  vgradp=ui(i,k)*psx(i)+vi(i,k)*psy(i)
675  os=os-vgradp*psi(i)*(dpmdps(i,k)-dpidps(i,k+1))- &
676  di(i,k)*(pm(i,k)-pi(i,k+1))
677  wi(i,k)=vgradp*psi(i)*dpmdps(i,k)+os
678  os=os-vgradp*psi(i)*(dpidps(i,k)-dpmdps(i,k))- &
679  di(i,k)*(pi(i,k)-pm(i,k))
680  enddo
681 !
682  enddo
683 !---
684  return
685  end subroutine
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:333
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:255
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:444