UPP v11.0.0
Loading...
Searching...
No Matches
GFSPOSTSIG.F
Go to the documentation of this file.
1
52subroutine 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! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225end 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
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:257
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:448
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:335