UPP (develop)
Loading...
Searching...
No Matches
GFSPOSTSIG.F
Go to the documentation of this file.
1
16!------------------------------------------------------------------------------------------
45subroutine rtsig(lusig,head,k1,k2,kgds,ijo,levs,ntrac,jcap,lnt2,me, &
46 h,p,px,py,t,u,v,d,trc,iret)
47
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
51 use omp_lib
52 implicit none
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
60!
61 integer idrt,io,jo,iidea
62! integer idrt,io,jo,mo3,mct,iidea,midea
63 integer(sigio_intkind):: irets
64! type(sigio_datm):: datm
65 type(sigio_dati) dati
66! type griddata
67! real,dimension(:,:),pointer :: datm
68! endtype griddata
69! type(griddata),dimension(:),pointer :: datatrc
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
74 integer idvm, klen
75 real pmean,sumq,xcp
76! integer, parameter :: latch=20
77! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78! Determine output grid
79 idrt = kgds(1)
80 if(kgds(1) == 0 .and. kgds(4) < 90000) idrt = 256
81 io = kgds(2)
82 jo = kgds(3)
83! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84! Read and transform surface fields
85 iret = 1
86 if (me == 0) then
87 print*,'Debug rtsig= ',lusig,k1,k2,ijo,kgds(1:20)
88 endif
89
90 idvm = head%idvm
91! jc = omp_get_num_threads()
92! write(*,*)' in RTSIG lnt2=',lnt2,' threads=',jc,' latch=',latch, &
93! ' jcap=',jcap,' io=',io,' jo=',jo,' ijo=',ijo
94!
95 if (k2 < k1) return
96
97 dati%i = 1 ! hs
98 dati%f => trisca(:,k1)
99 call sigio_rrdati(lusig,head,dati,irets)
100 if(irets /= 0) return
101
102! call sptez(0,jcap,idrt,io,jo,trisca(1,k1),h,1)
103! call sptez(0,jcap,idrt,io,jo,dats%hs,h,1)
104! call sptez(0,jcap,idrt,io,jo,dats%ps,p,1)
105! call sptezj(jcap,lnt2,1,idrt,io,jo,jc,trisca,h,latch,1)
106!
107 dati%i = 2 ! Surface pressure
108 dati%f => trisca(:,k1+1)
109 call sigio_rrdati(lusig,head,dati,irets)
110 if(irets /= 0) return
111!
112! call sptez(0,jcap,idrt,io,jo,trisca(1,k1),p,1)
113! call sptezj(jcap,lnt2,1,idrt,io,jo,jc,trisca,p,latch,1)
114!--
115 allocate(wrk(ijo,2))
116 call sptezm(0,jcap,idrt,io,jo,2,trisca(1,k1),wrk,1)
117 if( mod(idvm,10) < 2) then
118!$omp parallel do private(i)
119 do i=1,ijo
120 h(i) = wrk(i,1)
121 p(i) = 1.e3*exp(wrk(i,2))
122! p(i) = 1.e3*exp(p(i))
123 enddo
124 elseif(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) = 1000.*wrk(i,2)
129! p(i) = 1000.*p(i)
130 enddo
131 endif
132 if (allocated(wrk)) deallocate(wrk)
133
134 call sptezd(0,jcap,idrt,io,jo,trisca(1,k1+1),pmean,px,py,1)
135 iret = 0
136
137! if (k2 < k1) return
138
139! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140! Read and transform fields on levels k1 through k2
141 iret = 2
142 if (k2 >= k1) then
143 klen = k2-k1+1
144 do k=k1,k2
145 write(*,*)' retriving T for k=',k,' k1=',k1,' k2=',k2
146 dati%i = k + 2 ! Virtual Temperature or CpT
147 dati%f => trisca(:,k)
148
149 call sigio_rrdati(lusig,head,dati,iret)
150 enddo
151 call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),t(1,k1),1)
152! call sptezm(0,jcap,idrt,io,jo,klen,trisca,t,1)
153 do k=k1,k2
154 dati%i = levs + 2 + (k-1) * 2 + 1 ! Divergence
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 ! Vorticity
159 dati%f => triscb(:,k)
160 call sigio_rrdati(lusig,head,dati,irets)
161 if(irets /= 0) return
162 enddo
163 call sptezmv(0,jcap,idrt,io,jo,klen,trisca(1,k1),triscb(1,k1), &
164 u(1,k1),v(1,k1),1)
165 call sptezm(0,jcap,idrt,io,jo,klen,trisca(1,k1),d(1,k1),1)
166
167! call sptezm(0,jcap,idrt,io,jo,1,triscb,z(1,k),1)
168 write(*,*)' retriving d/z for k=',k,' k1=',k1,' k2=',k2
169! datm%z(3,:) = datm%z(3,:)+2*con_omega/sqrt(1.5)
170! call sptezm(0,jcap,idrt,io,jo,klen,datm%z,z,1)
171 write(*,*)' start get tracer'
172 do nt=1,ntrac
173 do k=k1,k2
174 dati%i = levs * (2+nt) + 2 + k ! Tracers starting with q
175 dati%f => trisca(:,k)
176 call sigio_rrdati(lusig,head,dati,irets)
177 enddo
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
180 enddo
181 !t=t/(1+con_fvirt*sh)
182 write(*,*)' end get tracer,idvm=',idvm,'ijo=',ijo,'ntrac=',ntrac
183!
184!-- get temp
185 if (mod(idvm/10,10) == 3) then ! Enthalpy case
186 allocate(cpi(0:ntrac))
187! write(*,*)'aft read sig, cpi=',head%cpi
188 cpi(0:ntrac) = head%cpi(1:ntrac+1)
189! write(*,*)'cpi=',cpi(0:ntrac)
190!$omp parallel do private(k,i,xcp,sumq,n)
191 do k=k1,k2
192 do i=1,ijo
193 xcp = 0.0
194 sumq = 0.0
195 do n=1,ntrac
196 if( cpi(n) /= 0.0 ) then
197 xcp = xcp + cpi(n)*trc(i,k,n)
198 sumq = sumq + trc(i,k,n)
199 endif
200 enddo
201 xcp = (1.-sumq)*cpi(0) + xcp
202 t(i,k) = t(i,k) / xcp ! Now g1 contains T
203 enddo
204 enddo
205 if (allocated(cpi)) deallocate(cpi)
206 else
207!$omp parallel do private(i,k)
208 do k=k1,k2
209 do i=1,ijo
210 t(i,k) = t(i,k) / (1+con_fvirt*trc(i,k,1)) !get temp from virtual temp
211 enddo
212 enddo
213 endif
214 endif
215! write(*,*)'end comput t'
216 iret=0
217! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
218end subroutine
219
220! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
248 subroutine modstuff(km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,&
249 pi,pm,om)
250 use sigio_module, only: sigio_modprd
251 implicit none
252 integer,intent(in):: km,idvc,idsl,nvcoord
253 real,intent(in):: vcoord(km+1,nvcoord)
254 real,intent(in):: ps,psx,psy
255 real,intent(in):: u(km),v(km),d(km)
256! real,intent(out):: pi(km+1),pm(km)
257 real*8, intent(out):: pi(km+1),pm(km)
258 real,intent(out):: om(km)
259 real*8 ps8,pm8(km),pd8(km),vcoord8(km+1,nvcoord)
260 real*8 dpmdps(km),dpddps(km),dpidps(km+1),pi8(km+1)
261 real vgradp,pd(km),px(km),py(km),os
262 integer k,iret,logk
263! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264 ps8=ps
265 vcoord8=vcoord
266 call sigio_modprd(1,1,km,nvcoord,idvc,idsl,vcoord8,iret,&
267 ps=(/ps8/),&
268 pm=pm8,pd=pd8,dpmdps=dpmdps,dpddps=dpddps)
269!
270!jw: has to be 8 real for wam
271 pi8(1)=ps
272 pm=pm8
273! pd=pd8
274 dpidps(1)=1.
275 do k=1,km
276 pi8(k+1)=pi8(k)-pd8(k)
277 dpidps(k+1)=dpidps(k)-dpddps(k)
278! if(pi(8)<0.) then
279! print *,'in modstuff,pi8=',pi8(k)
280! endif
281 enddo
282 pi=pi8
283!
284 os=0
285 do k=km,1,-1
286 vgradp=u(k)*psx+v(k)*psy
287 os=os-vgradp*ps*(dpmdps(k)-dpidps(k+1))-d(k)*(pm(k)-pi(k+1))
288 om(k)=vgradp*ps*dpmdps(k)+os
289 os=os-vgradp*ps*(dpidps(k)-dpmdps(k))-d(k)*(pi(k)-pm(k))
290 enddo
291 px=ps*dpmdps*psx
292 py=ps*dpmdps*psy
293 end subroutine
294
295!-------------------------------------------------------------------------------
327 subroutine modstuff2(im,ix,km,idvc,idsl,nvcoord,vcoord,ps,psx,psy,d,u,v,&
328 pi,pm,om,me)
329 use sigio_module, only : sigio_modprd
330 implicit none
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
338! real*8, allocatable :: ps8(:), pm8(:,:), pd8(:,:),dpmdps(:,:),dpddps(:,:), &
339! dpidps(:,:),pi8(:,:),vcoord8(:,:)
340! real, allocatable :: os(:)
341! real, allocatable :: pd(:,:),px(:,:), py(:,:), os(:)
342
343! real vgradpps
344
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)
348! real vgradpps,pd(im,km),px(im,km),py(im,km),os(im),tem
349 integer i,k,iret,logk
350! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351
352 ps8 = ps
353 vcoord8 = vcoord
354 call sigio_modprd(im,ix,km,nvcoord,idvc,idsl,vcoord8,iret, &
355 ps=ps8,pd=pd8,dpddps=dpddps,pm=pm8,dpmdps=dpmdps)
356
357!
358! if (me == 0) then
359! write(*,*)' pd8=',pd8(1,60:64)
360! write(*,*)' pm8=',pm8(1,60:64)
361! endif
362!jw: has to be 8 real for wam
363
364!$omp parallel do private(i)
365 do i=1,im
366 pi8(i,1) = ps(i)
367 dpidps(i,1) = 1.
368 os(i) = 0
369 pi(i,1) = pi8(i,1)
370 enddo
371 do k=1,km
372!$omp parallel do private(i)
373 do i=1,im
374 pi8(i,k+1) = pi8(i,k) - pd8(i,k)
375 dpidps(i,k+1) = dpidps(i,k) - dpddps(i,k)
376! if(pi(i,8)<0.) then
377! print *,'in modstuff,pi8=',pi8(i,k),' i=',i,' k=',k,' me=',me
378! endif
379 pi(i,k+1) = pi8(i,k+1)
380 pm(i,k) = pm8(i,k)
381 enddo
382 enddo
383!
384 do k=km,1,-1
385!$omp parallel do private(i,vgradpps)
386 do i=1,im
387 vgradpps = (u(i,k)*psx(i) + v(i,k)*psy(i)) * ps(i)
388
389 os(i) = os(i) - vgradpps*(dpmdps(i,k)-dpidps(i,k+1)) &
390 - d(i,k)*(pm(i,k)-pi(i,k+1))
391
392 om(i,k) = os(i) + vgradpps*dpmdps(i,k)
393
394 os(i) = os(i) - vgradpps*(dpidps(i,k)-dpmdps(i,k)) &
395 - d(i,k)*(pi(i,k)-pm(i,k))
396 enddo
397 enddo
398 end subroutine
399
400!-----------------------------------------------------------------------
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)
452 implicit none
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)
463 real tvcon,xcp,sumq
464 integer thermodyn_id,jn,js,is,in
465 integer jj,jjm,k,n,j,i,ij,lonb2
466! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
467! spectral transforms
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
471 lonb2=lonb*2
472 ij=lonb2*(j2-j1+1)
473 in=1
474 is=1+lonb
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)
485 if(j1==732)then
486 do k=1,km
487 do i=1,ijl
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)
492 end do
493 end do
494 end if
495 select case(mod(idvm,10))
496 case(0,1)
497 do i=1,ij
498 ps(i)=1.e3*exp(ps(i))
499 enddo
500 case(2)
501 do i=1,ij
502 ps(i)=1.e3*ps(i)
503 enddo
504 case default
505 do i=1,ij
506 ps(i)=1.e3*exp(ps(i))
507 enddo
508 end select
509! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
510 thermodyn_id=mod(idvm/10,10)
511 if (thermodyn_id == 3) then
512 do k=1,km
513 do i=1,ij
514 t(i,k) = t(i,k)/cpi(0) ! enthalpy (cpt/cpd)
515 end do
516 end do
517!
518 endif
519
520 call getomega(jcap,nc,km,idvc,idvm,idrt,idsl, &
521 nvcoord,vcoord,lonb,latb,ijl,j1,j2,1,sd,sps, &
522 ps,t,u,v,wi,pi,dpo)
523 if(j1==732)then
524 do k=1,km
525 do i=1,ijl
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)
529 end do
530 end do
531 end if
532 if(thermodyn_id /= 2)then
533! convert to surface pressure and temperature
534 if (thermodyn_id == 3) then
535 do k=1,km
536 do i=1,ij
537 xcp = 0.0
538 sumq = 0.0
539 do n=1,ntrac
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)
543 endif
544 enddo
545 t(i,k) = t(i,k)/((1.-sumq)*cpi(0)+xcp)
546 end do
547 end do
548
549 else
550 tvcon=(461.50/287.05-1.)
551 t(:,:) = t(:,:)/(1.+tvcon*q(:,1:km))
552 endif
553 end if
554 if(j1==732)then
555 do k=1,km
556 do i=1,ijl
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)
560 end do
561 end do
562 end if
563! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
564!----force tracers to be positive
565 if (chgq0 == 1) q = max(q, 0.0)
566! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
567! pass data to gfsdatao
568 do j=1,j2-j1+1
569 jn=j+j1-1
570 js=latb+1-jn
571 jn=(jn-1)*lonb
572 js=(js-1)*lonb
573 jj=j*lonb
574 jjm=(j-1)*lonb
575 do i=1,lonb
576 gfszs(i+jn) = zs(i+jjm)
577 gfsps(i+jn) = ps(i+jjm)
578 enddo
579 do i=1,lonb
580 gfszs(i+js) = zs(i+jj)
581 gfsps(i+js) = ps(i+jj)
582 enddo
583 do k=1,km
584 do i=1,lonb
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)
591 enddo
592 do i=1,lonb
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)
599 enddo
600 enddo
601 do k=1,km*ntrac
602 do i=1,lonb
603 gfsq(i+jn,k) = q(i+jjm,k)
604 enddo
605 do i=1,lonb
606 gfsq(i+js,k) = q(i+jj,k)
607 enddo
608 enddo
609 enddo
610 return
611 end
612!-----------------------------------------------------------------------
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)
642!!!!!
643 use sigio_module, only : sigio_modprd
644 implicit none
645
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)
652 real :: pi(ijn,km+1)
653 real :: os
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
658!----1. spectral transform
659 lonb2=lonb*2
660 ij=lonb2*(j2-j1+1)
661 in=1
662 is=1+lonb
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)
665
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)
668 psi8=psi
669 ti8=ti
670
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)
673 pm=pm8
674 pd=pd8
675
676 select case(mod(idvm,10))
677 case(0,1)
678 continue
679 case(2)
680 do i=1,ijn
681 psx(i)=psx(i)/(psi(i)*1.0e-3)
682 psy(i)=psy(i)/(psi(i)*1.0e-3)
683 enddo
684 case default
685 do i=1,ijn
686 psx(i)=psx(i)/psi(i)
687 psy(i)=psy(i)/psi(i)
688 enddo
689 end select
690
691!----3.omeda from modstuff
692 do i=1,ijn
693 pi(i,1)=psi(i)
694 dpidps(i,1)=1.
695 enddo
696 do k=1,km
697 do i=1,ijn
698 pi(i,k+1)=pi(i,k)-pd(i,k)
699 dpidps(i,k+1)=dpidps(i,k)-dpddps(i,k)
700 enddo
701 enddo
702 do i=1,ijn
703 os=0.
704 do k=km,1,-1
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))
711 enddo
712!
713 enddo
714!---
715 return
716 end subroutine
modstuff
subroutine modstuff(km, idvc, idsl, nvcoord, vcoord, ps, psx, psy, d, u, v, pi, pm, om)
Computes model coordinate dependent functions.
Definition GFSPOSTSIG.F:250
rtsig
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.
Definition GFSPOSTSIG.F:47
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)
Transforms sigma spectral fields to grid.
Definition GFSPOSTSIG.F:452
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)
Computes omega.
Definition GFSPOSTSIG.F:642
modstuff2
subroutine modstuff2(im, ix, km, idvc, idsl, nvcoord, vcoord, ps, psx, psy, d, u, v, pi, pm, om, me)
Computes model coordinate dependent functions.
Definition GFSPOSTSIG.F:329