UPP (upp-srw-2.2.0)
Loading...
Searching...
No Matches
GFSPOSTSIG.F
Go to the documentation of this file.
1
50
51! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94subroutine 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! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267end 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
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:299
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:96
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:501
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:691
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:378