subroutine wgradsupp(id7,NS,N,NL,iwmo,rla, rlo,p,z,t,td,u,v)
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
include 'tograds.inc'
dimension iwmo(N), rla(N), rlo(N), id7(7)
dimension p(N,NL),z(N,NL),t(N,NL),td(N,NL),u(N,NL),v(N,NL)
character *8 id
c UPPER-AIR grads file! no surface date!
c -------------------------------------------
if(NCTL.le.0) call getctls (-1, id7)
c -------------------------------------------
do 100 jfile=1, NCTL
print *,' file .... ',stn_ctl(jfile)
if (STN_CTL(jfile).ne.0) goto 101
100 continue
print *,' No ctl file !'
stop
101 continue
c ------------------------
id(8:8) = '\0'
iflag = 0
rt = 0
c ------------------------
do 200 k=1, NS
write ( id(1:5),'(i5)') iwmo(k)
write ( 91+jfile) id,rla(k),rlo(k),rt,NL,iflag
write ( 91+jfile)
+ (p(k,L),z(k,L),t(k,L),td(k,L),u(k,L),v(k,L),L=1,NL)
200 continue
rlat=0.
rlon = 0.
NLEV = 0
write ( 91+jfile) id,rlat,rlon,rt,NLEV,iflag
return
end
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
subroutine wgradsos(id7,rla, rlo,NNN,NPS,ips,itls,ils, VLS)
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
include 'tograds.inc'
c------------------------------------------------------------------
dimension id7(7),rla(NNN), RLO(NNN)
dimension IPS(NPS), ITLS(NPS), ILS(NPS), VLS(NNN,NPS)
dimension VHLP(1000*20)
character *8 CID
c**************************************************
if(NCTL.le.0) call getctls (-1, id7)
ihr = id7(5)
c**************************************************
c ----------------------
c ---------------------------------------------------
print *,' NNN :',ips, NNN, IPS,ITLS, ILS
do 50 jfile=1, NCTL
c ----------------------
if ( STN_CTL(jfile).eq.0) goto 50
c ----------------------
K1 = 1
do m = 1, JFILE - 1
NV = NVAR_CTL(jfile)
K1 = K1 + NV
enddo
NV = NVAR_CTL(jfile)
c ---------------------------------------------------
do 62 ip = 1, NPS
NPSF = 0
c ---------------------------------------------------
itl = itls(ip)
il = ils(ip)
ip1 = ips(ip)
do 600 K = K1, K1 + NV - 1 ! ctl variables
c print *,jfile, ip1, k, IP_CTL(k), ITL_CTL(k),IL_CTL(k)
if ( ip1.eq.IP_CTL(k).and.itl.eq.ITL_CTL(k).and.
+ il.eq.IL_CTL(k)) then
NPSF = 1
c print *,jfile, ' OK! :', k, ' p3:',ip1, itl,il
do i = 1, NNN
I1 = (i-1)*NV + K-K1+1
VHLP (I1) = VLS(i, ip)
enddo
goto 62
endif
600 continue ! ctl variables
62 continue ! income MODEL variable
c ---------------------------------------------------
IF ( NPSF.eq.0) goto 50
c ---------------------------------------------------
c writing!
rt = 0
do n = 1, NNN
cid (1:8)=' '
write (cid(1:5),'(i5)') N+3000
cid(8:8) = '\0'
write (91+jfile) cid, rla(n),rlo(n), rt, 1, 1 !IFLAG
M1 = (n-1)*NV + 1
write (*,*) n,NNN,cid, rla(n),rlo(n), (VHLP(m),m=m1, m1+NV-1)
write (91+jfile) (VHLP(m),m=M1,M1+NV-1 )
c write (*,*) NPSF, ' jfile=',jfile, (VHLP(m),m=M1,M1+NV-1 )
enddo
write (cid(1:5),'(i5)') NNN
write (91+jfile) cid, 0., 0., 0., 0,1
c write (*,*) cid, 0., 0., 0., 0,1
c ---------------------------------------------------
50 c o n t i n u e ! new ctl file
c ---------------------------------------------------
c*************************************************************
return
end
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c**************************************************
subroutine rgrads(id7,ip,itl,il,i9,j9,k9,zm,Vl1, NLRT)
c**************************************************
include 'tograds.inc'
dimension id7(7),zm(k9)
dimension vl1(i9,j9,k9)
c**************************************************
c print *,NCTL,id7,ip,itl,il,i9,j9,k9
c print *,zm
c**************************************************
if(NCTL.le.0) call getctls (0, id7)
c**************************************************
ihr = id7(5)
NLRT = 0
c print *,'RGRADS NCTL=',nctl,id7
c ----------------------
do 50 jfile=1, NCTL
c ----------------------
c ----------------------
if ( STN_CTL(jfile).eq.1) goto 50
c ----------------------
NB = index (ctln(jfile),' ') -1
c ----------------------
c print *,'RGRADS:',ctln(jfile)(1:NB),'xx'
c--------------------------------------------------
c access through the file name!
c------------------------------
if ( wctln(1:1).ne.' ') then
if ( wctln(1:NB).ne.ctln(jfile)(1:NB)) goto 50
endif
c--------------------------------------------------
RNX_CTL = NX_CTL(jfile)
RX0_CTL = X0_CTL(jfile)
RXI_CTL = XI_CTL(jfile)
RNY_CTL = NY_CTL(jfile)
RY0_CTL = Y0_CTL(jfile)
RYI_CTL = YI_CTL(jfile)
c print *,IHR, jfile, RNX_CTL, STN_CTL(jfile)
c ----------------------
NX = NX_CTL(jfile)
NY = NY_CTL(jfile)
NXNY=NX*NY
c if ( i9.ne.NX.or.j9.ne.NY) goto 50 ! call fmess ('check IM,JM!')
if ( i9*j9.ne.NX*NY) goto 50 ! call fmess ('check IM,JM!')
c ----------------------
c ----------------------
c print *,'RGRADS???:',IP,ITL,IL,zm,k9
call getrec
+ (jfile, id7, zm,k9,ip, 0, itl, il, NREC1,NREC2, NTIME)
c print *,' getrec:',nrec1,nrec2,NTIME
c ----------------------
if (NREC1.le.0) goto 50
c ----------------------
do 60 N=1, NREC1
c ----------------------
L = LMOD_CTL(N)
IR1 = IREC1_CTL(N)
read (91+JFILE,rec=IR1,err=60) (vall1 (m),m=1,NXNY)
c ----------------------
vmin = 100000
vmax = -vmin
ii = 0
c do 4 j=1, NY
c do 4 i=1, NX
do 4 j=1, J9
do 4 i=1, I9
ii = ii + 1
c print *,i,j,L,vall1(ii)
vmin = min(vmin,vall1(ii))
vmax = max(vmax,vall1(ii))
4 vl1 (i,j,L) = vall1 (ii)
c ----------------------
c IREC1_CTL(N) = 0
c IREC2_CTL(N) = 0
write(*,67) (id7(m),m=1,5),ir1,ctln(jfile)(NB-13+1:NB)
+ ,ip,itl,il, L
+, Vmin, Vmax
if ( vmin.lt.Vmax) NLRT = NLRT + 1
60 c o n t i n u e ! levels
50 c o n t i n u e ! file
67 format (' t:',5i2,'r:',i4,' f=',a13,4i5,' ',2f10.3)
return
end
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
subroutine wgrads(id7,ip1,ip2,iuv,itl,il,ar7,i9,j9,k9,zm,Vl1,Vl2)
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
include 'tograds.inc'
c include 'all.inc'
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c
c ARGUMENTs:
c------------------------------------------------------------------
c id7 - iyy, imm, idd, iss, ihr, imn, isc
c ip1 - parameter wmo-code number
c ip2 - parameter wmo-code for v comp. (otherwise 0)
c iuv - flag : 0 - parameter at "H" points
c 1 - u and v comp. but at "H" points
c 2 - u and v comp. at "V" points
c ITL - LEVEL-TYPE wmo-code number
c IL - LEVEL itself wmo-code number
c ar7 (area)...
c VAL1 [VAL2 ] - data
c------------------------------------------------------------------
dimension id7(7),ar7(7), zm(k9)
dimension vl1(i9,j9,k9), vl2(i9,j9,k9)
dimension vu(i9*j9), vv(i9*j9)
c**************************************************
if(NCTL.le.0) call getctls (1, id7)
ihr = id7(5)
c print *,' wgrads ... ctl NXNY=',j,NXNY, ' iuv=',iuv
c**************************************************
c ----------------------
do 50 jfile=1, NCTL
c ----------------------
if ( STN_CTL(jfile).eq.1) goto 50
c ----------------------
NXNY = NX_CTL(jfile) * NY_CTL(jfile)
NB = index (ctln(jfile),' ') -1
if ( iuv.eq.9) then
goto 55
IDXA = ar7(4)*1000
IDXACTL = XI_CTL(jfile)*1000
IDYA = ar7(7)*1000
IDYACTL = YI_CTL(jfile)*1000
if ( IDXA.ne.IDXACTL) goto 50
if ( IDYA.ne.IDYACTL) goto 50
if (ar7(2).ne.NX_CTL(jfile)) goto 50
if (ar7(3).ne.X0_CTL(jfile)) goto 50
if (ar7(5).ne.NY_CTL(jfile)) goto 50
if (ar7(6).ne.Y0_CTL(jfile)) goto 50
55 continue
endif
c ----------------------
call getrec
+ (jfile, id7, zm,k9,ip1, ip2, itl, il, NREC1,NREC2, NTIME)
if (NREC1.le.0) goto 50
c ----------------------
do 60 N=1, NREC1
L = LMOD_CTL(N)
IR1 = IREC1_CTL(N)
IR2 = IREC2_CTL(N)
if ( iuv.eq.9) then
print *,' direct write:',id7 ,ip1, itl, il,i9,k9
write (91+JFILE,rec=IR1) (((vl1 (i,j,k),i=1,i9),j=1,j9),k=1,k9)
goto 9898
endif
c call hcalmm (vl1(1,1,L),i9,j9,v1mn,v1mx)
c print *,' inp:',n,l, ir1, ir2, ' mn-mx:',v1mn, v1mx
if (ar7(1).eq.0.)
+call ll2ll (omn,omx,ar7,iuv,jfile,i9,j9,vl1(1,1,L),vl2(1,1,L))
if (ar7(1).eq.1.)
+call e1d2ll (omn,omx,ar7,iuv,jfile,i9,j9,vl1(1,1,L),vl2(1,1,L))
c if (ar7(1).eq.2.)
c +call e2d2ll (omn,omx,ar7,iuv,jfile,i9,j9,vl1(1,1,L),vl2(1,1,L))
if (ar7(1).eq.2.) then
call e221( iuv, vl1(1,1,L),vl2(1,1,L),i9,j9,vu,vv)
call e1d2ll (omn,omx,ar7,iuv,jfile,i9*j9,1,vu,vv)
endif
c if (ar7(1).eq.2.)
c +call e2d2ll (omn,omx,ar7,iuv,jfile,i9,j9,vl1(1,1,L),vl2(1,1,L))
c diagnostic!
7777 write (91+JFILE,rec=IR1) (vall1 (m),m=1,NXNY)
778 write(*,67) (id7(m),m=1,5),ir1,ctln(jfile)(NB-13+1:NB)
+ ,ip1,itl,il,zm(l)
+, omn,omx
67 format (' Wt=',5i2,'r:',i3,' ',a13,3i3,f6.0 ,2f6.0)
c -----------
9898 continue
IREC1_CTL(N) = 0
if (NREC2.le.0) goto 60
write(91+JFILE,rec=IR2) (vall2 (m),m=1,NXNY)
IREC2_CTL(N) = 0
60 c o n t i n u e ! levels
50 c o n t i n u e ! file
c*************************************************************
c*************************************************************
return
end
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
subroutine getctls (IRW, idat7)
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c IRW - READ-WRITE FLAG:
c IRW = 0 READ case date is untouched!
c = -1 Write-Observation
c = 1 Write-MODEL forecasts!
c character *10 d10
character *12 gdate
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
include 'tograds.inc'
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
dimension ihlp(20), rhlp(20), idat7(7)
character *150 str, dset
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
CALL findctl
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
call id72gdate (idat7, gdate)
print *,' grads:getctl!', gdate, ' NCTL=',NCTL
NPAR_CTL = 0
NLEV_CTL = 0
c---------------------------------------------------
do 400 k=1, NCTL
NOB = 0
NOB_CTL(k) = 0
STN_CTL(K) = 0
NVAR_CTL (k) = 0
c---------------------------------------------------
nb = index(CTLN(k),' ') - 1
if ( nb.lt.2) goto 400
c---------------------------------------------------
open (80,file=ctln(k)(1:nb),status='old')
open (81,file=ctln(k)(1:nb)//'.i')
JSTART = 0
c---------------------------------------------------
10 continue
c---------------------------------------------------
read (80,'(a)',end=99) str
if ( str(1:8).eq.'*obsdset') goto 14
if ( str(1:4).eq.'dset') goto 15
if ( str(1:7).eq.'*indset') goto 16
call s2lower(str)
c print *,str(1:60)
if ( index (str,'dtype stat').gt.0) STN_CTL(K) = 1
if ( str(1:4).eq.'tdef') goto 20
if ( str(1:4).eq.'xdef') goto 30
if ( str(1:4).eq.'ydef') goto 40
if ( str(1:4).eq.'zdef') goto 50
if ( str(1:4).eq.'vars') goto 60
if ( str(1:5).eq.'undef') goto 70
goto 10
c---------------------------------------------------
14 c o n t i n u e ! dset !
i1 = 10
m=i1
call mvwhile (' ','ne',m,+1,str)
i2 = m - 1
cVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
c +,NOB,XOB_CTL(NOBSMX),YOB_CTL(NOBSMX),ZOB_CTL(NOBSMX)
c +,JWMO_CTL(NOBSMX)
cVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
open (89,file=str(i1:i2),status='old',err=10)
141 read ( 89,*,end=142) JWMO, YOBS,XOBS,ZOBS
NOB = NOB + 1
JWMO_CTL(NOB) = JWMO
XOB_CTL(NOB) = XOBS
YOB_CTL(NOB) = YOBS
ZOB_CTL(NOB) = ZOBS
c print *,NOB,JWMO_CTL(NOB),YOB_CTL(NOB),XOB_CTL(NOB)
c . ,ZOB_CTL(NOB)
goto 141
142 close (89)
NOB_CTL(k) = NOB
print *,' NOB det',k,NOB_CTL(k)
c
goto 10
c---------------------------------------------------
15 c o n t i n u e ! dset !
i1 = 6
m=i1
call mvwhile (' ','ne',m,+1,str)
i2 = m - 1
lenfild = i2 - i1 +1
dset(1:lenfild) = str(i1:i2)
goto 10
c---------------------------------------------------
16 c o n t i n u e ! indset !
call clears (indset)
i1 = 9
m=i1
call mvwhile (' ','ne',m,+1,str)
i2 = m - 1
lenfil = i2 - i1 +1
indset(1:lenfil) = str(i1:i2)
print *,'indset=',indset(1:lenfil),'**'
goto 10
c------------------------------------
70 c o n t i n u e ! undef!
ips = 5
call s2r (str,ips, UNDEF_CTL(K))
goto 10
c---------------------------------------------------
20 c o n t i n u e ! time !
ips = 5
c------------------------------------
call s2i (str,ips, NT_CTL(K))
c------------------------------------
c IMPORTANT!
if ( ITF_CTL(k).ge.-1.or.IRW.eq.0) then
c------------------------------------
c get (respect) from the ctl file!
c------------------------------------
call gdate2hr (IHR,str,ips)
if (IHR.le.0) call fmess (ctln(K)//' TDEF')
c------------------------------------
else
c-----------------------------------------------------------
c get date from the program (model grb2ctl...)
c-----------------------------------------------------------
call datetohr
+ (idat7(3),idat7(2),idat7(1),idat7(4),0,IHR)
call s2i (str,ips, iDUM)
call s2i (str,ips, iDUM)
call s2i (str,ips, iDUM)
endif
c -----------------------------
c -----------------------------
c -----------------------------
c -----------------------------
IT0_CTL ( K ) = IHR
call s2i (str,ips, ITI_CTL(k))
c -----------------------------
if (index(str,'ye').gt.0) ITF1_CTL(k) = 6
if (index(str,'mo').gt.0) ITF1_CTL(k) = 5
if (index(str,'da').gt.0) ITF1_CTL(k) = 4
if (index(str,'hr').gt.0) ITF1_CTL(k) = 3
if (index(str,'mn').gt.0) ITF1_CTL(k) = 2
if (index(str,'sc').gt.0) ITF1_CTL(k) = 1
goto 10
c---------------------------------------------------
30 c o n t i n u e ! XDEF !
ips = 5
call s2i (str,ips, NX_CTL(K) )
call s2ar (str,ips, nh, rhlp)
X0_CTL(K) = rHLP(1)
XI_CTL(K) = rHLP(2)
c print *,k, ' xdef:',NX_CTL(k)
goto 10
c---------------------------------------------------
40 c o n t i n u e ! YDEF !
ips = 5
call s2i (str,ips, NY_CTL(K) )
call s2ar (str,ips, nh, rhlp)
c print *,k,' ydef:',NY_CTL(k)
Y0_CTL(K) = rHLP(1)
YI_CTL(K) = rHLP(2)
goto 10
c---------------------------------------------------
50 c o n t i n u e ! ZDEF !
ips = 4
call s2i (str,ips, NZ_CTL(K))
call s2ar (str,ips, NH0, rhlp)
do i =1, NH0
NLEV_CTL = NLEV_CTL + 1
NZE_CTL(k) = NLEV_CTL
Z_CTL(NLEV_CTL) = rhlp(i)
enddo
if (NH0.gt.NZ_CTL(k))
+ call fmess (ctln(k)(1:nb)//' ZDEF!')
if (NH0.eq.NZ_CTL(K)) goto 10
51 continue
read (80,'(a)',end=99) str
call s2lower(str)
ips = 1
call s2ar (str,ips, NH, rhlp)
do i = 1, NH
NLEV_CTL = NLEV_CTL + 1
NZE_CTL(k) = NLEV_CTL
if (NLEV_CTL.gt.NLEVMX) call fmess (ctln(k)(1:NB)//' ZDEF!')
Z_CTL(NLEV_CTL) = rhlp(i)
enddo
NH0 = NH0 + NH
if (NH0.gt.NZ_CTL(k)) call fmess (ctln(k)(1:NB)//' ZDEF!')
if ( NH0 .eq. NZ_CTL(K) ) goto 10
goto 51
c---------------------------------------------------
60 c o n t i n u e ! parameter !
c---------------------------------------------------
read (80,'(a)',end=99) str
c---------------------------------------------------
call s2lower(str)
if ( str(1:4).eq.'endv') goto 99
NPAR_CTL = NPAR_CTL + 1
ips = index (str,' ')
call s2i (str,ips, IPNL_CTL(NPAR_CTL))
call s2i (str,ips, IP_CTL (NPAR_CTL))
call s2i (str,ips, ITL_CTL (NPAR_CTL))
call s2i (str,ips, IL_CTL (NPAR_CTL))
IPF_CTL (NPAR_CTL) = K
IPREC_CTL(NPAR_CTL) = JSTART
c print *,NPAR_CTL,IPNL_CTL(NPAR_CTL),JSTART
JSTART = JSTART + max(1, IPNL_CTL (NPAR_CTL))
NREC_CTL(K) = JSTART
NVAR_CTL(k) = NVAR_CTL(k) + 1
c---------------------------------------------------
goto 60
c---------------------------------------------------
99 continue
c---------------------------------------------------
rewind (80)
80 continue
c print *,' REWIND80'
c update DATE!
read (80,'(a)',end=999) str
ips = 1
if ( str(1:4).eq.'tdef') then
if (IRW.ne.0.and.ITF_CTL(k).lt.-1) call srwd (str, ips, 4, gdate)
endif
m=150
call mvwhile (' ','eq',m,-1,str)
write (81,'(a)') str (1:m)//' '
goto 80
c---------------------------------------------------
999 continue
c---------------------------------------------------
close (81)
close (80)
c print *,' IRW=',IRW
call system ('mv '//ctln(k)(1:nb)//'.i '//ctln(k)(1:nb))
c open !
c if ( IRW.ge.0) then
if ( STN_CTL(K).ne.1) then
NXNY=NX_CTL(k)*NY_CTL(K)
print *,' NX*N,k,nob='
+,NXNY,k,NOB_CTL(k),' **',dset(1:lenfild),'**'
open ( 91+K,file=dset(1:lenfild),form='unformatted'
+ ,access='direct',recl=NXNY*1)
else
print *,' OBSERVATION ',dset(1:lenfild)
open ( 91+K,file=dset(1:lenfild),form='unformatted')
endif
400 continue
c print *,' grads:endgetctl!', NXNY
return
end
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
subroutine findctl
C NUOVO to ALOWWE various directories in USE!
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
include 'tograds.inc'
character *150 wplace,dirn, str
character *1 c1
call system ('pwd > wrk.grads')
do m=1,150
str(m:m)=' '
dirn(m:m)=' '
wplace(m:m)=' '
enddo
print *,' *** wrk?'
open (88,file='wrk.grads',status='old')
read (88,'(a)') wplace
nb = index (wplace,' ') -1
close (88)
print *,' ***',wplace(1:nb),'**'
open (88,file='activectl.LST',status='old',err=111)
c**********************************
10 continue
c**********************************
read (88,'(a)',end=99) str
c**********************************
i1 = index (str,'->')
caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
if ( i1.gt.0) then
if (index (str,wplace(1:nb)).eq.0) goto 10
c manage dirn!
i2 = 150
call mvwhile (' ','eq',i2,-1,str)
IDN = ichar(str(i2:i2))
if (IDN.ne.47) then
i2 = i2 + 1
str(i2:i2) = '/'
endif
i1 = i2
call mvwhile ('>','ne',i1,-1,str)
i1 = i1 + 1
k1 = 0
do k=i1,i2
c1 = str(k:k)
ic1 = ichar(c1)
if ( ic1.gt.32) then
k1 = k1 + 1
dirn(k1:k1) = c1
endif
enddo
kd = k1
endif
caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
ITIMEF =index (str,'=')
if ( ITIMEF .le.0) goto 10
if ( str(1:1).eq.'*') goto 10
c--------------------------------------
j1 = 1
call mvwhile (' ','eq',j1,+1,str)
j2 = j1
call mvwhile (' ','ne',j2,+1,str)
k1 = 0
do k=j1,j2
c1 = str(k:k)
ic1 = ichar(c1)
if ( ic1.gt.32) then
k1 = k1 + 1
str(k1:k1) = c1
endif
enddo
ks = k1
checkfile
open (89,file=dirn(1:kd)//str(1:ks),status='old',err=112)
close (89)
c print *,dirn(i1:i2)//str(1:nb)
ips = ITIMEF
c****************************
NCTL = NCTL + 1
call clears (ctln(NCTL))
ctln(NCTL)(1:kd) = dirn(1:kd)
ctln(NCTL)(kd+1:kd+ks) = str(1:ks)
ips = ips + 1
call s2i (str, ips, ITFLG)
ITF_CTL (NCTL) = ITFLG
print *,ITF_CTL(NCTL),' **',ctln(NCTL)(1:kd+ks),'**'
c----------------------------------------
goto 10
c----------------------------------------
99 continue
close (88)
return
111 print *,'STOP! NO current activectl.LST***'
stop
112 print *,'STOP:No file:',dirn(1:kd)//str(1:ks),'***'
stop
end
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
subroutine getrec
+ (NFILE,id7,zm,NZM,ip1,ip2,itl,il,NREC1,NREC2,NTP1)
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
include 'tograds.inc'
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
dimension ZM(NZM), id7(7)
NTP1 = 0
NREC1 = 0
NREC2 = 0
NIT = 1
if (ip2.gt.0) NIT = 2
do iter = 1, NIT
if (iter.eq.1) ip = ip1
if (iter.eq.2) ip = ip2
NREC= 0
NTIME=-1
do 200 k=1, NPAR_CTL
c------
c time!
c------
jfile = IPF_CTL(k)
if ( NFILE.ne.jfile) goto 200
c-------
c search!
c 0) time check!
c print *,' grads:getrec!', nzm, (zm(m),m=1,NZM),id7
call getNTIME (jfile,ID7,NTIME)
c print *,' grads:getrec!', NTIME
NTP1 = NTIME + 1
if (NTP1.lt.1.or.NTP1.gt.NT_CTL(jfile)) goto 200
PREC = NTIME * NREC_CTL (jfile) + IPREC_CTL(k)
c ------------------------------
if (IPNL_CTL(k).eq.0) goto 30 ! surface !
c ------------------------------
if(ip.ne.IP_CTL(k).or.itl.ne.ITL_CTL(k)) goto 200
c -------------------------------------------------
c find starting level!
LEV1 = 1
if (jfile.gt.1) LEV1 = NZE_CTL (jfile-1) + 1
LEV2 = NZE_CTL (jfile)
ILEV = min(LEV2 - LEV1, IPNL_CTL(k)-1)
do L = LEV1, LEV1 + ILEV
c print *,L,Z_CTL(L)
c ---------------------------
do 6 LMOD = 1, NZM
if ( ZM(LMOD).ne.Z_CTL (L)) goto 6
NREC = NREC + 1
LMOD_CTL(NREC) = LMOD
if(iter.eq.1) IREC1_CTL(NREC) = PREC+(L-LEV1+1)
if(iter.eq.2) IREC2_CTL(NREC) = PREC+(L-LEV1+1)
6 continue
5 enddo
goto 200
30 continue ! only surface!
c ---------------------------
c print *," getrec - surface "
if(ip.ne.IP_CTL(k).or.itl.ne.ITL_CTL(k).or.
+ il.ne.IL_CTL(k)) goto 200
c searc!
c-------
NREC = NREC + 1
LMOD_CTL(NREC) = 1
if(iter.eq.1) IREC1_CTL(NREC) = PREC + NREC
if(iter.eq.2) IREC2_CTL(NREC) = PREC + NREC
200 continue ! parameter
if(iter.eq.1) NREC1 = NREC
if(iter.eq.2) NREC2 = NREC
enddo ! iter
c print *,' grads:endgetrec!'
return
end
c----------------------------------------------------------------
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
subroutine getNTIME (jfile, id7, NTIME)
include 'tograds.inc'
dimension id7(7)
c----------------------------------------------------------------
NTIME = -999
if(ITI_CTL(JFILE).le.0) call fmess(ctln(JFILE)//' TDEF!')
call datetohr (id7(3),id7(2),id7(1),id7(4),id7(5),IHR)
if(ITF_CTL(jfile).ge.0) goto 50
c minutes and seconds!
IT0 = 0
if(ITF1_CTL(jfile).eq.6) then
IHLP=IHR - IT0_CTL(jfile)
ITI =365*24*ITI_CTL(jfile)
else if(ITF1_CTL(jfile).eq.5) then
IHLP=IHR - IT0_CTL(jfile)
ITI =28*24*ITI_CTL(jfile)
else if(ITF1_CTL(jfile).eq.4) then
IHLP=IHR - IT0_CTL(jfile)
ITI =24*ITI_CTL(jfile)
else if(ITF1_CTL(jfile).eq.3) then
IHLP=IHR - IT0_CTL(jfile)
ITI =ITI_CTL(jfile)
else if(ITF1_CTL(jfile).eq.2) then
IHLP=ID7(5)*60+ID7(6)
ITI =ITI_CTL(jfile)*60
else if(ITF1_CTL(jfile).eq.1) then
IHLP=ID7(5)*3600+ID7(6)*60 +id7(7)
ITI =ITI_CTL(jfile)*3600
else
call fmess (ctln(jfile)//' TDEF!')
endif
if( mod (IHLP,ITI).ne.0) return
NTIME = IHLP / ITI
return
50 continue
c ------------------------------------- ! speciale!
if ( ITF_CTL(jfile).ne.id7(5)) return
c ------------------------------------- ! speciale!
IDT = IHR - IT0_CTL (jfile)
NTIME = IDT/24
return
end
c----------------------------------------------------------------
subroutine ll2ll (vmin,vmax,ar7, iuv,jfile,i9,j9, val1, val2)
c----------------------------------------------------------------
include 'tograds.inc'
C-------------------------------------------------------------------
c dimension val1(imjm), val2(imjm)
dimension ar7(7),val1(i9,j9), val2(i9,j9)
common /first/ ifirst
c----------------------------------------------------------------
BLx(P,Q,X1,X2,X3,X4)=X1+P*(X2-X1)+Q*(X3-X1)+P*Q*(X1-X2-X3+X4)
c----------------------------------------------------------------
UND = UNDEF_CTL(jfile)
NPOI_CTL = 0
NX0 = ar7 (2)
WBD0 = ar7 (3)
DLMD0= ar7 (4)
NY0 = ar7 (5)
SBD0= ar7 (6)
DPHD0= ar7 (7)
ebd0 = wbd0 +(NX0-1)*DLMD0
rnbd0 = sbd0 +(NY0-1)*DPHD0
EBD0N = EBD0 + DLMD0
ICYCL = 0
if ( EBD0N.ge.360) EBD0N = EBD0N - 360.
if ( EBD0N.eq.WBD0) ICYCL = 1
NY = NY_CTL (jfile)
Y0 = Y0_CTL (jfile)
YI = YI_CTL (jfile)
NX = NX_CTL (jfile)
X0 = X0_CTL (jfile)
XI = XI_CTL (jfile)
defx = 1000*abs ( xi-dlmd0)
defy = 1000*abs ( yi-dphd0)
dex0 = abs(x0-WBD0)
dey0 = abs(Y0-SBD0)
vmin = 10000
vmax =-10000
ii = 0
c print *, XI, DLMD0, XI/DLMD0 , YI, DPHD0, YI/DPHD0
c print *, X0, WBD0, Y0, SBD0
c print *, NX, NX0, NY, NY0, defx, defy
if (NY.ne.NY0.or.NX.ne.NX0) goto 1000
if (dex0.ge.0001.or.dey0.ge.0001) goto 1000
if (defx.gt.1.0.or.defy.gt.1.0) goto 1000
c stop
ii = 0
do 5 j=1, NY
do 5 i=1, NX
ii = ii + 1
vall1 (ii) = val1(i,j)
vmin = min(vmin,vall1(ii))
vmax = max(vmax,vall1(ii))
if ( iuv.le.0) goto 5
vall2 (ii) = val2(i,j)
5 continue
NPOI_CTL = ii
if ( ifirst.le.0) ifirst = 0
ifirst = ifirst + 1
if(ifirst.eq.1) print *,' NO-GRADS interpolation! ',vmin,vmax
return
c ------------------------------------------------------
1000 continue
c stop
if ( ifirst.le.0) ifirst = 0
ifirst = ifirst + 1
if(ifirst.eq.1) print *,' YES-GRADS interpolation!',und
if(ifirst.eq.1) print *,' defx0-defy0:',defx, defy
if(ifirst.eq.1) print *,' dex0-dey0:',dex0, dey0
ii = 0
do 20 j= 1, NY
Y = Y0 + (j-1)*YI
do 55 i= 1, NX
X = X0 + (i-1)*XI
cpeja
CCNEWNEWNEW
if ( NOB_CTL(jfile).gt.0) then
if ( i.le.NOB_CTL(jfile)) then
X = XOB_CTL(i)
Y = YOB_CTL(i)
endif
endif
CCNEWNEWNEW
cpeja
ii = ii + 1
if (y.lt.sbd0) goto 55
if (y.eq.rnbd0) y=y-0.000001
if (y.gt.rnbd0) goto 55
if (x.lt.wbd0.or.x.gt.ebd0) then
x= x + 360
if (ICYCL.eq.0.and.(x.lt.wbd0.or.x.gt.ebd0)) goto 55
endif
NPOI_CTL = NPOI_CTL + 1
p = (x - wbd0)/dlmd0
q = (y - sbd0)/dphd0
i0 = p
j0 = q
p = p - i0
q = q - j0
i0 = i0 + 1
j0 = j0 + 1
i01 = i0 + 1
if ( ICYCL.eq.1.and.i01.eq.NX0+1) i01 = 1
j01 = j0 + 1
z1 = val1 (i0, j0)
z2 = val1 (i01,j0)
z3 = val1 (i01,j01)
z4 = val1 (i0 ,j01)
c vall1(ii)=blx(p,q, z1, z2, z3, z4)
vall1(ii)=blu(UND,p,q,z1,z2,z3,z4)
if(vall1(ii).ne.UND) vmin = min(vmin,vall1(ii))
if(vall1(ii).ne.UND) vmax = max(vmax,vall1(ii))
if ( iuv.le.0) goto 55
z1 = val2 (i0, j0)
z2 = val2 (i01,j0)
z3 = val2 (i01,j01)
z4 = val2 (i0, j01)
c vall2(ii)=blx(p,q, z1, z2, z3, z4)
vall2(ii)=blu (UND,p,q,z1, z2, z3, z4)
55 continue
20 continue
c stop
return
end
c----------------------------------------------------------------
subroutine llg2ll (vmin,vmax,ar7, iuv,jfile,i9,j9, val1, val2)
c----------------------------------------------------------------
include 'tograds.inc'
C-------------------------------------------------------------------
c dimension val1(imjm), val2(imjm)
dimension ar7(7),val1(i9,j9), val2(i9,j9)
common /gfirst/ NYG, AVDL, gaullat(400)
c----------------------------------------------------------------
BLx(P,Q,X1,X2,X3,X4)=X1+P*(X2-X1)+Q*(X3-X1)+P*Q*(X1-X2-X3+X4)
c----------------------------------------------------------------
UND = UNDEF_CTL(jfile)
vmin = 10000
vmax =-10000
NX0 = ar7 (2)
WBD0 = ar7 (3)
DLMD0= ar7 (4)
NY0 = ar7 (5)
SBD0= ar7 (6)
DPHD0= ar7 (7)
if ( NY0.ne.NYG) then
NYG = NY0
call gauss_lat_nmc
endif
ebd0 = wbd0 +(NX0-1)*DLMD0
rnbd0 = sbd0 +(NY0-1)*DPHD0
EBD0N = EBD0 + DLMD0
ICYCL = 0
if ( EBD0N.ge.360) EBD0N = EBD0N - 360.
c if ( EBD0N.gt.360) EBD0N = EBD0N - 360.
if ( EBD0N.eq.WBD0) ICYCL = 1
c print *,avdl, ' ICYCLE:', ICYCL,wbd0, ebd0, EBD0N
NY = NY_CTL (jfile)
Y0 = Y0_CTL (jfile)
YI = YI_CTL (jfile)
NX = NX_CTL (jfile)
X0 = X0_CTL (jfile)
XI = XI_CTL (jfile)
ii = 0
c***********************************************
do 70 J = 1, NY
c-----------------------------------------------
Y = Y0 + (j-1)*YI
JG = (Y - (-90.)) / avdl + 1
YG = gaullat (JG)
c -------------------
do while ( YG.gt.Y)
c -------------------
JG = JG - 1
YG = gaullat (JG)
c -------------------
enddo
c -------------------
DEY = gaullat(JG+1) - gaullat(JG)
q = (Y - YG)/DEY
c***********************************************
do 70 I = 1, NX
c***********************************************
ii = ii + 1
c-------------------------------------
X = X0 + (i-1)*XI
if (x.lt.wbd0.or.x.gt.ebd0) then
x= x + 360
c if (ICYCL.eq.0.and.(x.lt.wbd0.or.x.gt.ebd0)) goto 55
endif
p = (x - wbd0)/dlmd0
i0 = p
p = p - i0
c-------------------------------------
j0 = JG
j01 = j0+1
i0 = i0 + 1
i01 = i0 + 1
c----------------------------------------------
if ( ICYCL.eq.1.and.i01.eq.NX0+1) i01 = 1
c-----------------------------------------------
z1 = val1 (i0, j0)
z2 = val1 (i01,j0)
z3 = val1 (i01,j01)
z4 = val1 (i0 ,j01)
c vall1(ii)=blx(p,q, z1, z2, z3, z4)
vall1(ii)=blu(UND,p,q,z1,z2,z3,z4)
if(vall1(ii).ne.UND) vmin = min(vmin,vall1(ii))
if(vall1(ii).ne.UND) vmax = max(vmax,vall1(ii))
c if(ii.eq.NX*NY/2)print 55,x,y,p,q,i0,i01,j0,j01,vall1(ii)
70 continue
c----------------------------------------------------------------
55 format ( 2f9.4, 2f5.2, 4i4,f9.3)
c print *,' llg2ll:', vmin, vmax
c stop
return
end
c----------------------------------------------------------------
subroutine e1d2ll (vmin,vmax,ar7, iuv,jfile,i9,j9, val1, val2)
c----------------------------------------------------------------
include 'tograds.inc'
C-------------------------------------------------------------------
dimension ar7(7),val1(i9), val2(i9)
c----------------------------------------------------------------
BLx(P,Q,X1,X2,X3,X4)=X1+P*(X2-X1)+Q*(X3-X1)+P*Q*(X1-X2-X3+X4)
c----------------------------------------------------------------
UND = UNDEF_CTL(jfile)
TLM0D= ar7 (2)
WBD = ar7 (3)
DLMD = ar7 (4)
TPH0D= ar7 (5)
SBD = ar7 (6)
DPHD = ar7 (7)
IM= -wbd/dlmd +1.5
JM=2*(-sbd/dphd)+1.5
IMJM = im*jm - jm/2
KNE=IM
KNW=IM-1
NINC=2*IM-1
c----------------------------------------------------------------
dtr = .01745329
ctph0 = cos(tph0d*dtr)
stph0 = sin(tph0d*dtr)
c----------------------------------------------------------------
NY = NY_CTL (jfile)
Y0 = Y0_CTL (jfile)
YI = YI_CTL (jfile)
NX = NX_CTL (jfile)
X0 = X0_CTL (jfile)
XI = XI_CTL (jfile)
c print *,'e1d :',nx,x0,xi,ny,y0,yi, tlm0d
ii = 0
c------------------------------------>
vmin = 100000
vmax = -vmin
UND = UNDEF_CTL(jfile)
do 20 j= 1, NY
Y = Y0 + (j-1)*YI
do 55 i= 1, NX
X = X0 + (i-1)*XI
c--------------------------------------------
cpeja
CCNEWNEWNEWNEWNEW
if ( NOB_CTL(jfile).gt.0) then
if ( i.le.NOB_CTL(jfile)) then
X = XOB_CTL(i)
Y = YOB_CTL(i)
endif
endif
CCNEWNEWNEWNEWNEW
c--------------------------------------------
ii = ii + 1
vall1(ii)=UND
if( IUV.ne.0) vall2(ii)=UND
call ZTLLC(X,Y,TLM0D,DTR,CTPH0,STPH0,TLM,TPH)
if (tlm.le.wbd+2*dlmd ) goto 55
if (tlm.ge.-wbd-2*dlmd) goto 55
if (tph.le.sbd+2*dphd ) goto 55
if (tph.ge.-sbd-2*dphd) goto 55
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c tlm = tlm + .5 * mod ( im+1,2) * dlmd * dtr
c tph = tph + .5 * mod ( im+1,2) * dphd * dtr
call Zgcoef (ar7,tlm,tph, kh,ph,qh, kv,pv,qv)
if ( iuv.ne.2) then
KH1 = kh
KH2 = kh+ kne
KH3 = kh+ ninc
KH4 = kh+ knw
else
KH1 = kv
KH2 = kv + kne
KH3 = kv + ninc
KH4 = kv + knw
ph = pv
qh = qv
endif
c print *,ii,x,y,tlm,tph,ph,qh,kh1,kh2,kh3,kh4
if (iuv.gt.0) then
c v1=blx(ph,qh,val1(kh1),val1(kh2),val1(kh3),val1(kh4))
v1=blu(UND,ph,qh,val1(kh1),val1(kh2),val1(kh3),val1(kh4))
c v2=blx(ph,qh,val2(kh1),val2(kh2),val2(kh3),val2(kh4))
v2=blu(UND,ph,qh,val2(kh1),val2(kh2),val2(kh3),val2(kh4))
call ZRLTLW (X, Y, v1, v2, tlm0d,dtr,ctph0,stph0
+ , vall1(ii), vall2(ii))
if(vall1(ii).ne.UND) vmin = min(vmin,vall1(ii))
if(vall1(ii).ne.UND) vmax = max(vmax,vall1(ii))
else
c vall1(ii)=blx(ph,qh,val1(kh1),val1(kh2),val1(kh3),val1(kh4))
vall1(ii)=blu(UND,ph,qh,val1(kh1),val1(kh2),val1(kh3),val1(kh4))
if(vall1(ii).ne.UND) vmin = min(vmin,vall1(ii))
if(vall1(ii).ne.UND) vmax = max(vmax,vall1(ii))
endif
c print *,ii,x,y,tlm,tph,ph, qh,kh1,val1(kh1),vall1(ii)
55 continue
20 continue
c print *,'ende1dll', vmax, vmin
c stop
return
end
c----------------------------------------------------------------
subroutine e2d2ll (vmin,vmax,ar7, iuv,jfile,i9,j9, val1, val2)
c----------------------------------------------------------------
include 'tograds.inc'
c----------------------------------------------------------------
c include 'all.inc' ! !!!!!
C-------------------------------------------------------------------
c P A R A M E T E R
c & (IM=-WBD/DLMD+1.5,JM=-2*SBD/DPHD+1.5,imjm=im*jm-jm/2 )
C-------------------------------------------------------------------
c parameter (KNE=IM,KNW=IM-1, NINC=2*IM-1)
c----------------------------------------------------------------
c dimension val1(imjm), val2(imjm)
dimension ar7(7),val1(i9,j9), val2(i9,j9)
c----------------------------------------------------------------
BLx(P,Q,X1,X2,X3,X4)=X1+P*(X2-X1)+Q*(X3-X1)+P*Q*(X1-X2-X3+X4)
c----------------------------------------------------------------
TLM0D= ar7 (2)
WBD = ar7 (3)
DLMD = ar7 (4)
TPH0D= ar7 (5)
SBD = ar7 (6)
DPHD = ar7 (7)
IM= -(wbd/dlmd) +1.5
JM=2*(-sbd/dphd)+1.5
IMJM = im*jm - jm/2
KNE=IM
KNW=IM-1
NINC=2*IM-1
c----------------------------------------------------------------
dtr = .01745329
ctph0 = cos(tph0d*dtr)
stph0 = sin(tph0d*dtr)
c----------------------------------------------------------------
NY = NY_CTL (jfile)
Y0 = Y0_CTL (jfile)
YI = YI_CTL (jfile)
NX = NX_CTL (jfile)
X0 = X0_CTL (jfile)
XI = XI_CTL (jfile)
c print *,nx*ny,'e1d :',nx,x0,xi,ny,y0,yi, tlm0d
ii = 0
c------------------------------------>
vmin = 100000
vmax = -vmin
do 20 j= 1, NY
Y = Y0 + (j-1)*YI
do 55 i= 1, NX
X = X0 + (i-1)*XI
ii = ii + 1
c vall1 (ii) = undef
call ZTLLC(X,Y,TLM0D,DTR,CTPH0,STPH0,TLM,TPH)
if (tlm.le.wbd+2*dlmd ) goto 55
if (tlm.ge.-wbd-2*dlmd) goto 55
if (tph.le.sbd+2*dphd ) goto 55
if (tph.ge.-sbd-2*dphd) goto 55
call Zgcoef2 (ar7,tlm,tph, ih,jh,it,ph,qh, iv,jv,itv,pv,qv)
c tx= (tlm-wbd)/dlmd
c ty= (tph-sbd)/dphd
c call Zgcoef2P (tx,ty, ih,jh,it,ph,qh, iv,jv,itv,pv,qv)
c stop
m = it ! east
n = m -1 ! west
if ( iuv.eq.0) then
c print *,x,y,tx,ty,ih,itv,jh
vall1(ii)=blx
+(ph,qh,val1(ih,jh),val1(ih+m,jh+1),val1(ih+n,jh+1),val1(ih,jh+2))
else if ( iuv.eq.1) then
v1=blx
+(ph,qh,val1(ih,jh),val1(ih+m,jh+1),val1(ih+n,jh+1),val1(ih,jh+2))
v2=blx
+(ph,qh,val2(ih,jh),val2(ih+m,jh+1),val2(ih+n,jh+1),val2(ih,jh+2))
call ZRLTLW (X, Y, v1, v2, tlm0d,dtr,ctph0,stph0
+ , vall1(ii), vall2(ii))
else if ( iuv.eq.2) then
m = itv ! east
n = m -1 ! west
v1=blx
+(pv,qv,val1(iv,jv),val1(iv+n,jv+1),val1(iv+m,jv+1),val1(iv,jv+2))
v2=blx
+(pv,qv,val2(iv,jv),val2(iv+m,jv+1),val2(iv+n,jv+1),val2(iv,jv+2))
call ZRLTLW (X, Y, v1, v2, tlm0d,dtr,ctph0,stph0
+ , vall1(ii), vall2(ii))
endif
c print *,' *** XY ',x,y,tlm,tph,ih,jh,' val:',vall1(ii)
vmin = min(vmin,vall1(ii))
vmax = max(vmax,vall1(ii))
55 continue
20 continue
c print *,'ende2dll', vmax, vmin
c stop
return
end
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
SUBROUTINE ZRLTLW(ALMD,APHD,TPUS,TPVS,TLM0D,DTR,CTPH0,STPH0
& ,PUS,PVS)
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C
RELM=(ALMD-TLM0D)*DTR
SRLM=SIN(RELM)
CRLM=COS(RELM)
C
APH=APHD*DTR
SPH=SIN(APH)
CPH=COS(APH)
C
CC=CPH*CRLM
TPH=ASIN(CTPH0*SPH-STPH0*CC)
C
RCTPH=1./COS(TPH)
CRAY=STPH0*SRLM*RCTPH
DRAY=(CTPH0*CPH+STPH0*SPH*CRLM)*RCTPH
DC=DRAY*DRAY+CRAY*CRAY
PUS=(DRAY*TPUS+CRAY*TPVS)/DC
PVS=(DRAY*TPVS-CRAY*TPUS)/DC
C
RETURN
END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
subroutine Zgcoef2(ar7,tlm,tph,i2,j2,it,ph,qh,i2v,j2v,itv,pv,qv)
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c include 'all.inc'
c parameter
c +( im= -wbd/dlmd +1.5,jm=2*(-sbd/dphd)+1.5, imjm=im*jm-jm/2
c +, im1 = im-1, rim1=im-1)
dimension ar7(7)
c**********************************
WBD = ar7 ( 3)
DLMD= ar7 ( 4)
SBD = ar7 ( 6)
DPHD= ar7 ( 7)
im= -wbd/dlmd +1.5
jm=2*(-sbd/dphd)+1.5
im1 = im-1
rim1=im-1
k9 = 2*im - 1
c**********************************
X=(TLM-WBD)/DLMD
Y=(TPH-SBD)/DPHD
X1=0.5*( X+Y)
Y1=0.5*(-X+Y) ! +RIM1
I1=INT(X1)
J1=INT(Y1)
PH=X1-I1
QH=Y1-J1
JR=J1 ! -IM1
I2=I1-JR
J2=I1+JR
i2 = (i2-1)/2+1
it = mod(j2+1,2)
c KH=J2*IM-J2/2+(I2+2)/2
C-----
XV1=X1-.5
YV1=Y1-.5
I1=INT(XV1)
J1=INT(YV1)
PV=XV1-I1
QV=YV1-J1
JR=J1 ! -IM1
I2V=I1-JR
J2V=I1+JR
i2v = (i2v-1)/2+1
itv = mod(j2v+1,2)
c KV=(J2+1)*IM-(J2+1)/2+(I2+1)/2
C-----
return
end
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
subroutine Zgcoef (ar7,tlm,tph, kh,ph,qh, kv,pv,qv)
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c include 'all.inc'
c parameter
c +( im= -wbd/dlmd +1.5,jm=2*(-sbd/dphd)+1.5, imjm=im*jm-jm/2
c +, im1 = im-1, rim1=im-1)
dimension ar7(7)
c**********************************
WBD = ar7 ( 3)
DLMD= ar7 ( 4)
SBD = ar7 ( 6)
DPHD= ar7 ( 7)
im= -wbd/dlmd +1.5
jm=2*(-sbd/dphd)+1.5
im1 = im-1
rim1=im-1
c**********************************
X=(TLM-WBD)/DLMD
Y=(TPH-SBD)/DPHD
X1=0.5*( X+Y)
Y1=0.5*(-X+Y)+RIM1
I1=INT(X1)
J1=INT(Y1)
PH=X1-I1
QH=Y1-J1
JR=J1-IM1
I2=I1-JR
J2=I1+JR
KH=J2*IM-J2/2+(I2+2)/2
C-----
XV1=X1-.5
YV1=Y1-.5
I1=INT(XV1)
J1=INT(YV1)
PV=XV1-I1
QV=YV1-J1
JR=J1-IM1
I2=I1-JR
J2=I1+JR
KV=(J2+1)*IM-(J2+1)/2+(I2+1)/2
C-----
return
end
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
SUBROUTINE ZTLLC(ALMD,APHD,TLM0D,DTR,CTPH0,STPH0,TLM,TPH)
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
RELM=(ALMD-TLM0D)*DTR
SRLM=SIN(RELM)
CRLM=COS(RELM)
C
APH=APHD*DTR
SPH=SIN(APH)
CPH=COS(APH)
C
CC=CPH*CRLM
ANUM=CPH*SRLM
DENOM=CTPH0*CC+STPH0*SPH
C-----------------------------------------------------------------------
TLM=ATAN2(ANUM,DENOM)/DTR
TPH=ASIN(CTPH0*SPH-STPH0*CC)/DTR
C-----------------------------------------------------------------------
RETURN
END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c********************************************
c string library!
c********************************************
c subroutine s2i ( str, ips, I4)
c subroutine s2r ( str, ips, R4)
c subroutine s2iass ( str, ips, sstr, i4)
c subroutine s2rass ( str, ips, sstr, r4)
c********************************************
c subroutine s2ai ( str, ips, NI4, I4A)
c subroutine s2ar ( str, ips, NI4, R4A)
c********************************************
c subroutine s2s1b ( str, str1, lstr1)
c subroutine s21b ( str, lstr1)
c********************************************
c subroutine s2sb2z ( str, str1)
c subroutine s2b2z ( str)
c********************************************
c subroutine s2slower ( str, str1)
c subroutine s2lower ( str)
c********************************************
c subroutine sPs ( str, ips, str1,ips1)
c********************************************
c subroutine fmess ( str )
c********************************************
c
c********************************************
subroutine mindex (k,idir,sstr,lstr,str)
c********************************************
character (*) str
character (*) sstr
c lstr = len (str)
lsstr = len (sstr)
i1 = k
i2 =lstr-lsstr+1
is = 1
if ( idir.lt.0) then
i1 = lstr-lsstr+1
i2 = k
is = -1
endif
c print *,' LSTR=',lstr,lsstr,' MINDEX:',i1,i2,is
print *,' mindex:',i1,i2,is
do k=i1,i2,is
if ( index(str(k:k+lsstr-1), sstr).gt.0) goto 7
enddo
k=0
7 return
end
c********************************************
subroutine sindex (k,idir,sstr,lstr,str)
c********************************************
character (*) str
character (*) sstr
c lstr = len (str)
lsstr = len (sstr)
i1 = k
i2 = LSTR -lsstr+1! lstr-lsstr+1
is = 1
if ( idir.lt.0) then
i1 = k -lsstr+1
i2 = 1
is = -1
endif
print *,'**',sstr(1:lsstr),'**',lsstr,' sX:',i1,i2,is
do k=i1,i2,is
if ( index(str(k:k+lsstr-1), sstr).gt.0) goto 7
enddo
k=0
7 return
end
c********************************************
subroutine mvwhile (c1,OPT,m,idirec, str)
c********************************************
character (*) str
character *1 c1
character *2 OPT
lstr = len (str)
if ( m.le.0) call fmess (str)
if ( m.gt.lstr) call fmess (str)
c print *,'mvwhie ',OPT,m,lstr,idirec
c print *,'**',str(1:lstr),'**'
c if ( idirec.lt.0) then
c idirec = -1
c else
c idirec = 1
c endif
ic1 = ichar(c1)
is = ichar(str(m:m))
c print *,'mvwhie ',OPT,ic1, is
if ( OPT.eq.'ne') then
do while (is.ne.ic1)
m = m + idirec
c if ( m.le.0) call fmess (str)
if ( m.le.0.or.m.gt.LSTR) return
is = ichar(str(m:m))
enddo
else if ( OPT.eq.'le') then
do while (is.le.ic1)
m = m + idirec
c if ( m.le.0) call fmess (str)
if ( m.le.0) return
is = ichar(str(m:m))
enddo
else if ( OPT.eq.'lt') then
do while (is.lt.ic1)
m = m + idirec
c if ( m.le.0) call fmess (str)
if ( m.le.0) return
is = ichar(str(m:m))
enddo
else if ( OPT.eq.'ge') then
do while (is.ge.ic1)
m = m + idirec
c if ( m.le.0) call fmess (str)
if ( m.le.0) return
is = ichar(str(m:m))
enddo
else if ( OPT.eq.'gt') then
do while (is.gt.ic1)
m = m + idirec
c if ( m.le.0) call fmess (str)
if ( m.le.0) return
is = ichar(str(m:m))
enddo
else if ( OPT.eq.'eq') then
c print *,' -------eq:',m,is
do while (is.eq.ic1)
m = m + idirec
c if ( m.le.0) call fmess (str)
if ( m.le.0) return
is = ichar(str(m:m))
c print *,m,is
enddo
endif
return
end
c********************************************
c---------------------------------------------------------
subroutine s2i ( str, ips, I4)
c---------------------------------------------------------
character (*) str
character *1 c1
lstr = len (str)
i = ips
in = 0
ie = 0
10 c1 = str (i:i)
if ( c1.ge.'0'.and.c1.le.'9'.and.in.eq.0) in = i
if ((c1.lt.'0'.or .c1.gt.'9').and.in.ne.0) goto 30
if ( i.eq.lstr) goto 20
i = i + 1
goto 10
20 ips = -1
return
30 ips = i
c read ( str(in:i-1), '(i8)') i4
c print *,'s2i:',str(in:i-1),'**'
read ( str(in:i-1), *) i4
if(str(in-1:in-1).eq.'-' ) i4 = -i4
if ( ips.eq.lstr) ips = 0
return
end
c---------------------------------------------------------
subroutine s2r ( str, ips, R4)
c---------------------------------------------------------
character (*) str
character *1 c1
lstr = len (str)
i = ips
in = 0
ie = 0
10 c1 = str (i:i)
if ( c1.ge.'0'.and.c1.le.'9'.and.in.eq.0) in = i
if ( c1.eq.'.' .and.in.eq.0) in = i
if ( (c1.eq.' ').and.in.ne.0) goto 30
if ( (c1.eq.',').and.in.ne.0) goto 30
if ( (c1.eq.';').and.in.ne.0) goto 30
if ( (c1.eq.':').and.in.ne.0) goto 30
c if ( (c1.eq.'.').and.in.ne.0.and.in.ne.i) goto 30
if ( i.eq.lstr) goto 20
i = i + 1
goto 10
20 ips = -1
return
c SOLUTION!
30 ips = i
c read ( str(in:i-1), '(f12.4)') R4
read ( str(in:i-1), *) R4
if(str(in-1:in-1).eq.'-' ) r4 = -r4
c print *,in,i-1,' **',str(in:i-1)
if ( ips.eq.lstr) ips = 0
return
end
c---------------------------------------------------------
subroutine s2iass ( str, ips, sstr, i4)
c---------------------------------------------------------
character (*) str
character (*) sstr
ips = -1
lstr = len (str)
lsstr = len (sstr)
i = index ( str, sstr )
if ( i.eq.0) return
ips = i + lsstr - 1
call s2i ( str, ips, i4)
return
end
c---------------------------------------------------------
subroutine s2rass ( str, ips, sstr, r4)
c---------------------------------------------------------
character (*) str
character (*) sstr
ips = -1
lstr = len (str)
lsstr = len (sstr)
i = index ( str, sstr )
if ( i.eq.0) return
ips = i + lsstr - 1
call s2r ( str, ips, r4)
return
end
c---------------------------------------------------------
subroutine s2ai ( str, ips, NI4, I4A)
c---------------------------------------------------------
character (*) str
dimension I4A (1)
lstr = len (str)
NI4 = 0
call s2i ( str, ips, i4)
do while ( ips.gt.0)
NI4 = NI4 + 1
I4A (NI4) = i4
call s2i ( str, ips, i4)
enddo
return
end
c---------------------------------------------------------
subroutine s2ar ( str, ips, NI4, R4A)
c---------------------------------------------------------
character (*) str
dimension R4A (1)
c dimension R4A (100)
c dimension (*) R4A
lstr = len (str)
NI4 = 0
call s2r ( str, ips, r4)
do while ( ips.gt.0)
NI4 = NI4 + 1
R4A (NI4) = r4
call s2r ( str, ips, r4)
enddo
return
end
c---------------------------------------------------------
subroutine s2s1b ( str, str1, lstr1)
c---------------------------------------------------------
character (*) str
character (*) str1
lstr = len (str)
c leading blank!
k0 = 1
do while (str(k0:k0) .le. ' ')
k0 = k0 + 1
enddo
k1 = 0
do 6 k = k0, lstr
c if ( str(k:k+1).eq.' ') goto 6
if ( str(k:k).le.' '.and.str(k+1:k+1).le.' ') goto 6
k1 = k1 + 1
str1 (k1:k1) = str(k:k)
6 continue
lstr1 = k1
return
end
c---------------------------------------------------------
subroutine s21b ( str, lstr1)
c---------------------------------------------------------
character (*) str
lstr = len (str)
c leading blank!
k0 = 1
ic1=ichar(str(k0:k0))
do while (ic1.le.32)
k0 = k0 + 1
ic1=ichar(str(k0:k0))
enddo
k1 = 0
do 6 k = k0, lstr
c if ( str(k:k+1).eq.' ') goto 6
ic1=ichar(str(k:k))
ic2=ichar(str(k+1:k+1))
c if ( str(k:k).le.' '.and.str(k+1:k+1).le.' ') goto 6
if ( ic1.le.32.and.ic2.le.32) goto 6
k1 = k1 + 1
str (k1:k1) = str(k:k)
6 continue
lstr1 = k1
return
end
c---------------------------------------------------------
subroutine s2sb2z ( str, str1)
c---------------------------------------------------------
character (*) str
character (*) str1
lstr = len (str)
c leading blank!
do k =1, lstr
str1 (k:k) = str(k:k)
if ( str(k:k).eq.' ') str1(k:k) = '0'
enddo
return
end
c---------------------------------------------------------
subroutine s2b2z ( str)
c---------------------------------------------------------
character (*) str
lstr = len (str)
c leading blank!
do k =1, lstr
if ( str(k:k).eq.' ') str(k:k) = '0'
enddo
return
end
c---------------------------------------------------------
subroutine s2slower ( str, str1)
c---------------------------------------------------------
character (*) str
character (*) str1
lstr = len (str)
do k =1, lstr
str1(k:k) = str(k:k)
ic1 = ichar (str (k:k))
if ( ic1.ge.65.and.ic1.le.97) str1(k:k)=char(ic1+32)
enddo
return
end
c---------------------------------------------------------
subroutine s2lower ( str)
c---------------------------------------------------------
character (*) str
lstr = len (str)
do k =1, lstr
ic1 = ichar (str (k:k))
if ( ic1.ge.65.and.ic1.le.90) str(k:k)=char(ic1+32)
enddo
return
end
c---------------------------------------------------------
subroutine scps ( str, ips, str1)
c---------------------------------------------------------
character (*) str
character (*) str1
l1 = len (str1)
str (ips:ips + l1 - 1) = str1 (1:l1)
ips = ips + l1 - 1
return
end
c---------------------------------------------------------
subroutine sPs ( str, ips, str1,ips1)
c---------------------------------------------------------
character (*) str
character (*) str1
if ( ips.le.0) return
l0 = len (str )
if ( ips + ips1 .gt. l0 ) call fmess
+ ('Length string problem:')
do k = 1, ips1
k0 = k + ips
str (k0:k0) = str1(k:k)
enddo
ips = k0
return
end
c---------------------------------------------------------
c subroutine sUpdate ( str, ips, NW, sstr)
subroutine sRwd ( str, ips, NW, sstr)
c---------------------------------------------------------
c replase NW word in str by sstr!
c--------------------------------
character *1500 str1
character (*) str
character (*) sstr
call clears (str1)
call s2s1b (str,str1, l1)
c call s21b (sstr,lss)
lss = len (sstr)
NWC = 0
do k=1, l1
if ( str1(k:k).eq.' ') then
NWC = NWC + 1
if ( NWC.eq.NW-1) k1 = k - 1
if ( NWC.eq.NW.or.k.eq.l1) then
k2 = k + 1
goto 6
endif
endif
enddo
c---------------------------------------------------------
call fmess ('FE:word number is wrong'//str1(1:l1))
c---------------------------------------------------------
6 continue
c print *,k1,k2, '****',str1(k1:k2), '** NW=',nw
call clears (str)
c----------------
c the first part!
c----------------
ips = 1
call scps (str, ips, str1(1:k1))
ips = ips + 2
call scps (str, ips, sstr(1:lss))
ips = ips + 2
call scps (str, ips, str1(k2:l1))
return
end
c---------------------------------------------------------
subroutine clears ( str)
c---------------------------------------------------------
character (*) str
lstr = len (str)
do k =1, lstr
str (k:k) = ' '
enddo
return
end
c---------------------------------------------------------
subroutine ra2sY (str,ips,NDEC,NSPc,NI4,RA4)
c---------------------------------------------------------
character (*) str
character * 15 hlp
dimension RA4 (1)
c NINROW = 65/(NSPc+1) - 1
NINROW = 6
ips1 = ips
c call clears ( str )
c ) if ndec = 0 ** no dot in string!!!
CCC ips = 1
c ips = 15
do 20 l = 1, NI4
hlp = ' '
write ( hlp(1:15),* ) RA4(l)
idot = index ( hlp(1:15),'.')
ke = idot + NDEC
k = 1
do while (hlp(k:k).eq.' ')
k = k + 1
enddo
lenh = ke - k + 1
call s2b2z ( hlp (k:ke))
str (ips:ips+lenh-1) = hlp (k:ke)
c ips = ips + lenh + 1
last = ips+lenh + 0
c------------------------------------------
if ( mod(l,NINROW).eq.0) then
str(last:last) = char(10)
c------------------------------------------
ips = last + ips1 + 1
else
c------------------------------------------
ips = ips + NSPc ! space for number!
endif
c------------------------------------------
20 continue
return
end
c---------------------------------------------------------
c---------------------------------------------------------
subroutine ra2s (str, ips, NDEC, NI4, RA4)
c---------------------------------------------------------
character (*) str
character * 15 hlp
dimension RA4 (1)
c call clears ( str )
c ) if ndec = 0 ** no dot in string!!!
CCC ips = 1
c ips = 1
do 20 l = 1, NI4
hlp = ' '
write ( hlp(1:15),* ) RA4(l)
idot = index ( hlp(1:15),'.')
ke = idot + NDEC
k = 1
do while (hlp(k:k).eq.' ')
k = k + 1
enddo
lenh = ke - k + 1
call s2b2z ( hlp (k:ke))
str (ips:ips+lenh-1) = hlp (k:ke)
ips = ips + lenh + 1
20 continue
return
end
c---------------------------------------------------------
subroutine m3tom2 ( str, ips, im2)
c---------------------------------------------------------
character (*) str
character *36 cmon
data cmon /'janfebmaraprmayjunjulaugsepoctnovdec'/
i = index (cmon, str(ips:ips+2))
if ( i.le.0) call fmess ('FE:MONTH PROBLEM')
im2 = (i-1)/3 + 1
ips = ips + 3
return
end
c---------------------------------------------------------
c---------------------------------------------------------
c subroutine dateGrtoHR ( str, ips, T0)
c---------------------------------------------------------
c call s2i ( str
c return
c end
c---------------------------------------------------------
c---------------------------------------------------------
subroutine fmess ( str )
c---------------------------------------------------------
character (*) str
l1 = len (str)
print *,'GRADS ERROR:',str(1:l1)
stop
end
c ******************************************************
subroutine hrtodate
+ ( ihr0, iss,ihh, idd, imm, iyy,ihhf)
c ******************************************************
c This routine returns from passed hours in current
c century - proper forecast time according to
c started and retrieve forecast time!
c Pejanovicz (Apr 97.)
c ******************************************************
dimension mont (12)
dimension ipr (21)
data ipr /464592,499656,534720,569784,604848,639912,674976
. ,710040,745104,780168,815232,850296,885360 ! 2000
. ,920424,955488,990552,1025616,1060680,1095744
c----------------------------------------------------
. ,1130808,1165872/
cod 1952-2032
data mont /31,28,31,30,31,30,31,31,30,31,30,31/
mont(2) = 28
ihr = ihr0
iprest=0
do k=1,21
if ( ihr.ge.ipr(k).and.ihr.lt.ipr(k)+24) then
iprest=1
iyy=1952+(k-1)*4
imm=12
idd=31
ihhf = ihr-ipr(k)
endif
enddo
if(iprest.eq.1) then
return
endif
c print *," date tohr ... ihr=",ihr
if ( ihr.ge.920424.and.ihr.lt.920448) then
iyy=2004
imm=12
idd=31
ihhf = ihr-920424
c print *," now interv ",ihr ! dodato za 31.12.2004!!!
return
endif
ihrp = ihr - iss - ihh
ipd = ihrp/24
ipyfg = ipd/365
subs = 1
jtr = 0
12 rpyfg = ipd/365.
c print *," ipyfg =",ipyfg, rpyfg,subs
nyy29 = (ipyfg-1)/4
ryy29 = (rpyfg-subs)/4.0
nyy29 = ryy29
c print *,ipd," nny29 ", nyy29,ryy29
if (ipd.ge.366) nyy29 = nyy29 + 1
nyy28 = (ipd - nyy29*366 ) / 365
c-------------------------
iyy = nyy28 + nyy29
c-------------------------
ipdcy = ipd - nyy29*366 - nyy28*365
ipdcm = 0
if ( mod(iyy,4).eq.0) mont(2) = 29
do j=0, 11
if ( ipdcm+mont(j+1).gt.ipdcy) goto 60
ipdcm = ipdcm + mont(j+1)
enddo
60 continue
c-------------------------
imm = j + 1
idd = ipdcy - ipdcm + 1
c if ( mod(iyy,4).eq.0.and.imm.eq.12 ) subs = 1.5
if ( mod(iyy,4).eq.0.and.imm.ge.11 ) subs = 1.5
jtr = jtr +1
if ( jtr.eq.1) goto 12
c-------------------------
ihhf = ihr - nyy29*366*24 - nyy28*365*24 - ipdcy*24 -iss
c print '(5i3)',iyy,imm,idd,iss,ihhf
c print *,' iyy=',iyy,' imm=',imm,' idd=',idd,iss,ihhf
C2000.
c if ( mont(2).eq.29.and.imm.eq.12.and.idd.gt.6) IHR = IHR + 24
iyy =iyy + 1900
return
end
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c PREPROCESING !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c ******************************************************
subroutine datetohr( idd, imm, iyy1,iss,ihh, ihr)
c ******************************************************
c This routine calculates passed hours in current 20-th
c century where iss is started time and ihh required
c forecast time!
c Pejanovicz (Apr 97.)
c ******************************************************
dimension mont (12)
data mont /31,28,31,30,31,30,31,31,30,31,30,31/
mont(2) = 28
c new lines (29aug98)
iyy = mod (iyy1,100)
if ( iyy.lt.50) iyy = iyy + 100
c--------------------
nyy29 = (iyy-1)/4
if ( iyy.gt.0) nyy29 = nyy29+1
nyy28 = iyy - nyy29
ipd = nyy29*366+nyy28*365 + idd - 1
if ( mod(iyy,4).eq.0) mont(2) = 29
do j = 1,imm -1
ipd = ipd + mont(j)
enddo
c if ( imm.gt.2.and.mod(iyy,4).eq.0) ipdcy = ipdcy + 1
ihr = ipd * 24 + iss + ihh
c print*, ' ihr = ',ihr
c print '(5i3,i8)',iyy,imm,idd,iss,ihh, ihr
return
end
c000000000000000000000000000000000000000000000000000000000000
subroutine cdate2id7 (cdate8,id7)
character *8 cdate8
dimension id7 (7)
id7 (5) = 0
id7 (6) = 0
id7 (7) = 0
do j = 1, 4
j1 = (j-1)*2 + 1
read (cdate8(j1:j1+1),'(i2)') id7 (j)
enddo
return
end
c000000000000000000000000000000000000000000000000000000000000
subroutine gdate2hr (JHR, str,ips)
c000000000000000000000000000000000000000000000000000000000000
character (*) str
character *36 cmon
data cmon /'janfebmaraprmayjunjulaugsepoctnovdec'/
JHR = -1
call mvwhile (' ','le',ips,+1,str)
c do while (str(ips:ips).le.' ')
c ips = ips + 1
c enddo
call s2i (str,ips,ISS)
call s2i (str,ips,IDD)
idx =index ( cmon,str(ips:ips+2))
if (idx.eq.0) return
imm =(idx-1) / 3 + 1
call s2i (str,ips,IYY)
call datetohr (idd,imm,iyy,iss,00, JHR)
return
end
c000000000000000000000000000000000000000000000000000000000000
subroutine id72gdate(idat7,gdate)
c000000000000000000000000000000000000000000000000000000000000
dimension idat7(7)
character *12 gdate
character *36 cmon
data cmon /'janfebmaraprmayjunjulaugsepoctnovdec'/
write (gdate (1:2), '(i2)'), idat7(4)
gdate (3:3) = 'Z'
write (gdate (4:5), '(i2)'), idat7(3)
iy = idat7(1)
C2000.
if (iy.gt.99) goto 1999
if (iy.lt.50) then
iy = iy + 2000
else
iy = iy + 1900
endif
1999 continue
c print *,'wwwwwwww:',idat7
c print *,'wwwwwwww:iy=',iy
im = idat7(2)
i = (im-1)*3 + 1
gdate (6:8) = cmon (i:i+2)
write (gdate(9:12),'(i4)' ) iy
call s2b2z (gdate)
return
end
c---------------------------------------------
subroutine hcalmm ( X,i1,j1,vmin,vmax)
dimension X(i1,j1)
vmin = 100000
vmax = -vmin
do j=1, j1
do i=1, i1
vmin = min (vmin, X(i,j))
vmax = max (vmax, X(i,j))
enddo
enddo
return
end
c---------------------------------------------
subroutine e221( iuv,val1,val2,im,jm,vu,vv)
c---------------------------------------------
dimension val1(im,jm),val2(im,jm),vu(im*jm),vv(im*jm)
c flag iuv = 0 ! scalar filed stored in the H-points
c flag iuv = 1 ! vector fields stored in the H-points
c flag iuv = 2 ! vector fields stored in the H-points
c --------------------
if (iuv.eq.2) goto 2000
c --------------------
k=0
do j =1, jm
do i =1, im - mod (j+1,2)
k = k + 1
vu (k) = val1(i,j)
if ( iuv.eq.1) vv (k) = val2(i,j)
enddo
enddo
return
c do j=1,jm
c do i=1,im
c if(.not.(mod(j,2).eq.0.and.i.eq.im))then
c k=k+1
c vu (k) = val1(i,j)
c if ( iuv.eq.1) vv (k) = val2(i,j)
c endif
c enddo
c enddo
c-----------------------------------------
2000 continue
c-----------------------------------------
k=0
do j =1, jm
do i =1, im - mod (j,2)
k = k + 1
vu (k)=val1(I,J)
vv (k)=val2(I,J)
enddo
enddo
c paran broj tacaka!
c do j=1,jm
c do i=1,im
c if(.not.(mod(j,2).eq.1.and.i.eq.im))then
c k=k+1
c vu (k)=val1(I,J)
c vv (k)=val2(I,J)
c endif
c enddo
c enddo
end
c------------------------------------------
function blu (UNDEF,p,q,V1,V2,V3,V4)
c------------------------------------------
SS = 0.
VV = 0.
if ( V1.ne.UNDEF) then
S = (1.-p)*(1.-q)
SS = SS + S
VV = VV + V1* S
endif
if ( V2.ne.UNDEF) then
S = p*(1.-q)
SS = SS + S
VV = VV + V2* S
endif
if ( V3.ne.UNDEF) then
S = p*q
SS = SS + S
VV = VV + V3* S
endif
if ( V4.ne.UNDEF) then
S = (1.-p)*q
SS = SS + S
VV = VV + V4* S
endif
if (SS.ne.0.) then
blu = VV /SS
else
blu = UNDEF
endif
c print *,SS,VV,blu,v1,v2,v3,v4
return
end
c------------------------------------------
c------------------------------------------
subroutine Zgcoef2P ( X,Y, ih,jh,it,p,q, iv,jv,itv,pv,qv)
c------------------------------------------
PI=3.1415926
a = PI/180.*45.
s2 = sqrt (2.)
rs2 = 1./s2
c print *,' Enter x-y'
c read *,x,y
xp=rs2 * ( x*cos(a) + y*sin(a))
yp=rs2 * (-x*sin(a) + y*cos(a))
ip =int (xp)
if (xp.lt.0) ip = ip -1
p = xp - ip
jp =int (yp)
if (yp.lt.0) jp = jp -1
q = yp - jp
c print *,ip,jp,xp,yp,p,q
xp = float(ip)
yp = float(jp)
x0=s2 * ( xp*cos(-a) + yp*sin(-a))
y0=s2 * (-xp*sin(-a) + yp*cos(-a))
ih = nint(x0)+1
ih = (ih-1)/2+1
jh = nint(y0)+1
it = mod(jh,2) - 1
c print *,ih,jh,it,p,q
c print *,nint(x0),nint(y0),p,q
iv = ih
jv = jh
itv = it
pv = p
qv = q
return
end
c SUBROUTINE readpMW3FT33x(NSFLAG,OUT)
SUBROUTINE MW3FT33x(IRET, NSFLAG,OUT)
C SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: W3FT33 THICKEN THINNED WAFS GRIB GRID 37-44
C PRGMMR: RALPH PETTERSON ORG: W/NMCXX DATE: 94-11-13
C
C ABSTRACT: SUBROUTINE THICKENS ONE THINNED WAFS GRIB GRID TO A
C REAL ARRAY OF 5329 NUMBERS (73,73) 1.25 DEGREE GRID.
C
C PROGRAM HISTORY LOG:
C 94-??-?? RALPH PETERSON
C 94-11-07 R.E.JONES ADD DOC BLOCK, CHANGE CALL TO 3
C PARAMETERS. REPLACE COS WITH TABLE
C LOOKUP.
C
C USAGE: CALL W3FT33(AIN, OUT, NSFLAG)
C INPUT ARGUMENT LIST:
C AIN - REAL 3447 WORD ARRAY WITH UNPACKED THINNED WAFS
C GRIB TYPE 37-44.
C NSFLAG - INTEGER = 1 AIN IS WAFS GRIB GRID 37-40 N. HEMI.
C = -1 AIN IS WAFS GRIB GRID 41-44 S. HEMI.
C
C OUTPUT ARGUMENT LIST:
C OUT - REAL (73,73) WORD ARRAY WITH THICKENED WAFS GRIB
C GRID 37-44.
C
C REMARKS: THE POLE POINT FOR U AND V WIND COMPONENTS WILL HAVE ONLY
C ONE POINT. IF YOU NEED THE POLE ROW CORRECTED SEE PAGE 9 SECTION
C 1 IN OFFICE NOTE 388. YOU NEED BOTH U AND V TO MAKE THE
C CORRECTION.
C
C ATTRIBUTES:
C LANGUAGE: SiliconGraphics 5.2 FORTRAN 77
C MACHINE: SiliconGraphics IRIS-4D/25, 35, INDIGO, Indy
C
C$$$
C
PARAMETER (NX=73,NY=73)
PARAMETER (NIN=3447)
C
REAL AIN(NIN)
REAL OUT(NX,NY)
C
INTEGER IPOINT(NX)
C
c SAVE
cKGDS = 0 65535 73 0 -120000 128 90000 -30000 65535 1250 64 0 0 0 0 0 0 0 0 33
c73 73 73 73 73 73 73 73 72 72 72 71 71 71 70
c70 69 69 68 67 67 66 65 65 64 63 62 61 60 60
c59 58 57 56 55 54 52 51 50 49 48 47 45 44 43
c42 40 39 38 36 35 33 32 30 29 28 26 25 23 22
c20 19 17 16 14 12 11 9 8 6 5 3 2
DATA IPOINT/
& 73, 73, 73, 73, 73, 73, 73, 73, 72, 72, 72, 71, 71, 71, 70,
& 70, 69, 69, 68, 67, 67, 66, 65, 65, 64, 63, 62, 61, 60, 60,
c & 70, 69, 69, 68, 68, 67, 66, 65, 65, 64, 63, 62, 61, 60, 60,
& 59, 58, 57, 56, 55, 54, 52, 51, 50, 49, 48, 47, 45, 44, 43,
& 42, 40, 39, 38, 36, 35, 33, 32, 30, 29, 28, 26, 25, 23, 22,
& 20, 19, 17, 16, 14, 12, 11, 9, 8, 6, 5, 3, 2/
C
iret = -1
open (77,file='dump',form='unformatted',status='old',err=900)
read (77 ,err=901,end=901) AIN
close (77)
c print *,' READOK!'
NXM = NX - 1
FNXM = FLOAT(NXM)
SPCG = 90.0 / FNXM
IE = 0
RADS = 3.14159/180.
DO J = 1,NY
C ALAT = FLOAT(J-1) * SPCG
C IF (NSFLAG.LT.0) ALAT=90.-ALAT
C NPOINT = IFIX(2.0+(90./SPCG)*COS(ALAT*RADS))
IF (NSFLAG.LT.0) THEN
NPOINT = IPOINT(74-J)
ELSE
NPOINT = IPOINT(J)
END IF
C IF (NPOINT.GT.NX) NPOINT = NX
IS = IE + 1
IE = IS + NPOINT - 1
C IF (IE.GT.NIN) STOP
C WRITE (*,*) ALAT,NPOINT,IS,IE
DPTS = (FLOAT(NPOINT)-1.) / FNXM
PW = 1.0
PE = PW + DPTS
OUT(1,J) = AIN(IS)
OUT(NX,J) = AIN(IE)
VALW = AIN(IS)
VALE = AIN(IS+1)
DVAL = (VALE-VALW)
DO I = 2,NXM
WGHT = PE -FLOAT(IFIX(PE))
OUT(I,J) = VALW + WGHT * DVAL
PW = PE
PE = PE + DPTS
IF (IFIX(PW).NE.IFIX(PE)) THEN
IS = IS + 1
VALW = VALE
VALE = AIN(IS+1)
DVAL = (VALE - VALW)
END IF
ENDDO
ENDDO
iret = 0
900 continue
RETURN
901 continue
close (77)
RETURN
END
c------------------------------------------------------------
subroutine readdump ( IRET,IYREV,IM,JM, VAL)
dimension VAL (IM,JM)
c------------------------------------------------------------
IRET = -1
open (77,file='dump',form='unformatted',status='old',err=901)
if ( IYREV.eq.1) then
read (77,err=900,end=900) ((VAL(i,j),i=1,IM),J=JM,1,-1)
else
read (77,err=900,end=900) ((VAL(i,j),i=1,IM),J=1,JM)
endif
close (77)
IRET = 0
return
900 close (77)
901 return
end
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c GrADS-INTERFACE Library;
c By Pejanovicz Aug-98 Malta
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
subroutine bsslz1(bes,n)
c
c
implicit real*8(a-h,o-z)
dimension bes(n)
dimension bz(50)
c
data pi/3.14159265358979d0/
data bz / 2.4048255577d0, 5.5200781103d0,
$ 8.6537279129d0,11.7915344391d0,14.9309177086d0,18.0710639679d0,
$ 21.2116366299d0,24.3524715308d0,27.4934791320d0,30.6346064684d0,
$ 33.7758202136d0,36.9170983537d0,40.0584257646d0,43.1997917132d0,
$ 46.3411883717d0,49.4826098974d0,52.6240518411d0,55.7655107550d0,
$ 58.9069839261d0,62.0484691902d0,65.1899648002d0,68.3314693299d0,
$ 71.4729816036d0,74.6145006437d0,77.7560256304d0,80.8975558711d0,
$ 84.0390907769d0,87.1806298436d0,90.3221726372d0,93.4637187819d0,
$ 96.6052679510d0,99.7468198587d0,102.888374254d0,106.029930916d0,
$ 109.171489649d0,112.313050280d0,115.454612653d0,118.596176630d0,
$ 121.737742088d0,124.879308913d0,128.020877005d0,131.162446275d0,
$ 134.304016638d0,137.445588020d0,140.587160352d0,143.728733573d0,
$ 146.870307625d0,150.011882457d0,153.153458019d0,156.295034268d0/
nn=n
if(n.le.50) go to 12
bes(50)=bz(50)
do 5 j=51,n
5 bes(j)=bes(j-1)+pi
nn=49
12 do 15 j=1,nn
15 bes(j)=bz(j)
return
end
subroutine gauss_lat_nmc
c
implicit real*8(a-h,o-z)
dimension a(500)
c-------------------------------------------------
real*4 gaullat(400) , avdl
common /gfirst/ NYG, avdl, gaullat
c-------------------------------------------------
k = NYG
c
c save
c
c k = 190
esp=1.d-14
c=(1.d0-(2.d0/3.14159265358979d0)**2)*0.25d0
fk=k
kk=k/2
call bsslz1(a,kk)
do 30 is=1,kk
xz=cos(a(is)/sqrt((fk+0.5d0)**2+c))
iter=0
10 pkm2=1.d0
pkm1=xz
iter=iter+1
if(iter.gt.10) go to 70
do 20 n=2,k
fn=n
pk=((2.d0*fn-1.d0)*xz*pkm1-(fn-1.d0)*pkm2)/fn
pkm2=pkm1
20 pkm1=pk
pkm1=pkm2
pkmrk=(fk*(pkm1-xz*pk))/(1.d0-xz**2)
sp=pk/pkmrk
xz=xz-sp
avsp=abs(sp)
if(avsp.gt.esp) go to 10
a(is)=xz
30 continue
if(k.eq.kk*2) go to 50
a(kk+1)=0.d0
pk=2.d0/fk**2
do 40 n=2,k,2
fn=n
40 pk=pk*fn**2/(fn-1.d0)**2
50 continue
do 60 n=1,kk
l=k+1-n
a(l)=-a(n)
60 continue
c
radi=180./(4.*atan(1.))
gstart = 0.
dga = 0.
NN = 0.
do 211 n=1,k
gaullat(n)=acos(a(n))*radi-90.0
dg = gaullat(n) - gstart
if ( n.gt.1) then
dga = dga + dg
NN = NN + 1
endif
c print '(i5,f9.3, f9.5)',n, gaullat(n),dg
gstart = gaullat(n)
211 continue
degla = dga/NN
c print *,'gaussian lat (deg) for jmax=',degla
c44 print *,' UNESI LAT:'
c read *,rglat
c k = (rglat - (-90.)) / degla + 1
c print '(4f9.3)',gaullat(k-1),gaullat(k),rglat,gaullat(k+1)
c print *,(gaullat(n),n=1,k)
c
avdl = degla
return
70 write(6,6000)
6000 format(//5x,14herror in gauaw//)
stop
end
c---------------------------------------------------------
subroutine s2capital( str)
c---------------------------------------------------------
character (*) str
lstr = len (str)
do k =1, lstr
ic1 = ichar (str (k:k))
c if ( ic1.ge.65.and.ic1.le.90) str(k:k)=char(ic1+32)
if ( ic1.ge.97.and.ic1.le.122) str(k:k)=char(ic1-32)
enddo
return
end
c---------------------------------------------------------
subroutine fconvert (ID,JD,IS,IC,sstr, string,LSTR)
c---------------------------------------------------------
character (*) sstr
character (*) string
character *1 z1,z2
character *5 c5,c6
z1 = char(39)
z2 = char(34)
call clears(string)
call clears(c5)
call clears(c6)
L = 0
ie = 13
ie = 28
string(1:ie) = "/usr/X11R6/bin/convert -fill"
c collor 1 yellow , 2 red ............................
if ( IC.eq.1) then
n = 6
ie = ie +1 + n
string(ie+1-n:ie) = "yellow"
else if ( IC.eq.2) then
n = 3
ie = ie +1 + n
string(ie+1-n:ie) = "red"
else if ( IC.eq.3) then
n = 5
ie = ie +1 + n
string(ie+1-n:ie) = "green"
else if ( IC.eq.4) then
n = 4
ie = ie +1 + n
string(ie+1-n:ie) = "blue"
else if ( IC.eq.5) then
n = 5
ie = ie +1 + n
string(ie+1-n:ie) = "black"
else
n = 5
ie = ie +1 + n
string(ie+1-n:ie) = "white"
endif
c pointsize
n = 10
ie = ie +1 + n
string(ie+1-n:ie) = "-pointsize"
n = 3
ie = ie+1 +n
write(string(ie+1-n:ie),'(i3)') IS
c draw
n = 11
ie = ie +1 + n
string(ie+1-n:ie) = "-draw "//z1//"text"
write(c5,'(i5)') ID
call s21b ( c5, L5)
print *,ID,'c5=',c5(1:L5),'**'
write(c6,'(i5)') JD
call s21b ( c6, L6)
print *,JD,'c6=',c6(1:L6),'**'
n = L5+1+L6
ie = ie+1 + n
string(ie+1-n:ie) = c5(1:L5)//","//c6(1:L6)
c text
lss= index(sstr,' ') -1
n = lss+1+1+1
ie = ie+1 +n
string(ie+1-n:ie) = z2//sstr(1:lss)//z2//z1
c print *,' sstr=',sstr(1:LSS),'*'
call s21b ( string, L)
LSTR = L
print *,'*',string(1:L),'*'
return
end
c---------------------------------------------------------
subroutine fcompose (I,J,IV,sstr, string,L)
c call fcompose (IS,JS,JVR,"sim$$$.gif")
c string1(1:20)="composite -geometry "
c string2(1:20)="composite -geometry "
c---------------------------------------------------------
character (*) sstr
character (*) string
character *1 z1,z2
character *5 c5,c6
z1 = char(39)
z2 = char(34)
call clears(string)
call clears(c5)
call clears(c6)
L = 0
c geomtery
c=======================================================
Cpeja ie = 19
Cpeja string(1:ie) = "composite -geometry"
ie = 58
string(1:ie) =
& "/usr/local/ImageMagick-5.4.4/utilities/composite -geometry"
C 1 2 3 4 5
C 123456789012345678901234567890123456789012345678901234567890
write(c5,'(i5)') I
call s21b ( c5, L5)
write(c6,'(i5)') J
call s21b ( c6, L6)
n = 1+L5+1+L6
ie = ie+1 + n
string(ie+1-n:ie) = "+"//c5(1:L5)//"+"//c6(1:L6)
c=======================================================
c text
LS0= index(sstr,' ') -1
jd = index(sstr,'$$$')
write ( sstr(jd:jd+2),'(i3)') IV
LS = 0
do k=1, LS0
if ( ichar(sstr(k:k)).eq.32) then
else
LS = LS +1
sstr(LS:LS) = sstr(k:k)
endif
enddo
n = LS
ie = ie+1 +n
string(ie+1-n:ie) = sstr(1:LS)
c print *,' sstr=',sstr(1:LSS),'*'
call s21b ( string, L)
return
end
c000000000000000000000000000000000000000000000000000000000000
subroutine dayinweek ( id7, jday )
c000000000000000000000000000000000000000000000000000000000000
dimension id7 (7)
data IHR0 /898368/
data iday0 /3/
call datetohr( id7(3), id7(2), id7(1),id7(4),id7(5), ihr)
nd = (IHR-IHR0)/24.
if ( mod(IHR-IHR0,24).lt.0) then
jday=mod(iday0-1+nd,7)
else
jday = mod ( iday0+nd,7 )
endif
if ( jday.lt.0) jday = 7+jday
return
end
c=============================================
subroutine rsynop( str,LS,iwmo,id7,data)
dimension id7(7),data(19)
c=============================================
character *(LS) str
character *3 c3
LST= len ( str)
id7(5) = 0
id7(6) = 0
id7(7) = 0
c3 = ","//char(92)//"N"
i = 1
call s2i (str,i, iwmo)
call s2i (str,i, id7(1))
call s2i (str,i, id7(2))
if ( id7(2).le.0) id7(2) = -id7(2)
call s2i (str,i, id7(3))
if ( id7(3).le.0) id7(3) = -id7(3)
call s2i (str,i, id7(4))
call s2i (str,i, idum)
call s2i (str,i, idum)
i = i + 1
do 211 k=1, 19
data(k) = -99.
if ( i.ge.LST) goto 211
if ( index(str(i:i+2),c3) .le.0) then
call s2r (str,i, data(k))
else
i = i + 3
endif
211 continue
return
end
c=============================================
subroutine wsynop( str,LS,iwmo,id7,data)
dimension id7(7),data(19)
c=============================================
character *(LS) str
character *3 c3
character *2 c2
c3 = char(32)//char(92)//"N"
c LS = len (str)
id7(5) = 0
id7(7) = 0
id7(7) = 0
i = 1
c====================================
write ( str(1:5),'(i5.5)') iwmo
c====================================
IYY = id7(1)
if ( IYY.lt.40) then
IYY = 2000+IYY
else
if( IYY.lt.100) IYY = 1900+IYY
endif
str(6:6) = char(44) ! ,
str(7:7) = char(34) ! "
str(8:26) = "1997-02-13+00:00:00"
write ( str(8:11),"(i4.4)") IYY
write ( str(13:14),"(i2.2)") id7(2)
write ( str(16:17),"(i2.2)") id7(3)
write ( str(19:20),"(i2.2)") id7(4)
str(27:27) = char(34) ! "
c str(28:28) = char(44) ! "
i = 29
do k=1, 19
rv = data(k)
iv = data(k)
takeiv = 1
if ( k.ge.4.and.k.le.9) takeiv = 0
if ( iv.eq.-99) then
str(i:i+2) = c3
c i = i + 3
i = i + 5
else
c str(i:i) = char(44)
if (takeiv.eq.1) write(str(i:i+4),"(i5)") iv
if (takeiv.eq.0) write(str(i:i+5),"(f6.1)")rv
i = i + 8
endif
enddo
LS = len (str)
call s21b (str,LS)
do k=1, LS
if ( ichar(str(k:k)).eq.32) str(k:k) = char(44)
if ( str(k:k).eq."+") str(k:k) = " "
enddo
if (str(LS:LS).eq."," )LS = LS -1
return
end
c===========================================================
subroutine pisigs(LFN,rlat,rlon,ISM,SZ,str)
c===========================================================
c
character (*) str
character *1 c1
character *4 csz
character *50 ctitle
L1 = len(str)
if ( L1.gt.5) L1 = index(str," ")
c1 = char(34)
write(LFN,*) '"query ll2xy ',rlon,rlat,'"'
write(LFN,*) 'x=subwrd(result,1)'
write(LFN,*) 'y=subwrd(result,2)'
if ( ISM.eq.0) then
ctitle=
. c1//"draw string "
. //c1//"x"//c1//" "//c1//"y"//c1//" "//str(1:L1)//" "//c1
write(LFN,"(a)")ctitle
c write(*,"(a)")ctitle
else if ( ISM.gt.0) then
c ctitle = draw mark 1 "x" "y" 0.08"'
write ( csz,"(f4.2)") SZ
ctitle=
. c1//"draw mark "//char(48+ISM)//" "
. //c1//"x"//c1//" "//c1//"y"//c1//" "//csz//" "//c1
write(LFN,"(a)")ctitle
endif
return
end
subroutine jday(IYY,IMM,IDD,julday)
dimension monthday(11)
DATA monthday/31,28,31,30,31,30,31,31,30,31,30/
if (mod(iyy,4).eq.0)monthday(2)=29
julday=IDD
do m=1,imm-1
julday=julday+monthday(m)
enddo
return
end
subroutine date2jday(IYY,IMM,IDD,julday)
dimension monthday(12)
DATA monthday/31,28,31,30,31,30,31,31,30,31,30,31/
if (mod(iyy,4).eq.0)monthday(2)=29
if ( idd.lt.1.or.idd.gt.monthday(imm)) then
julday = 0
return
endif
julday=IDD
do m=1,imm-1
julday=julday+monthday(m)
enddo
return
end
subroutine jday2date(IYY,IMM,IDD,julday)
c this routine calculates IMM (month), IDD, day in month
c input :IYY (year), julday (julian day)
c output :IMM (month in year IYY)
c :IDD (day in month)
dimension monthday(11)
DATA monthday/31,28,31,30,31,30,31,31,30,31,30/
if (mod(iyy,4).eq.0)monthday(2)=29
i1= 1
do m=1,11
i2 = i1 + monthday(m)-1
if ( julday.ge.i1.and.julday.le.i2) then
goto 5
endif
c print *,m,i1,i2
i1 = i2+1
enddo
5 continue
IMM = m
IDD = julday - i1 + 1
return
end
subroutine hbox(ar7,IM,JM,rlatk,rlonk,ih,jh)
PARAMETER (PI=3.141592654, s2d=1.4142136)
PARAMETER (DTR=PI/180.0)
dimension ar7(7)
tlm0d = AR7(2)
wbd = AR7(3)
dlmd = AR7(4)
tph0d = AR7(5)
sbd = AR7(6)
dphd = AR7(7)
TPH0=TPH0D*DTR
CTPH0=COS(TPH0)
STPH0=SIN(TPH0)
almd = rlonk
aphd = rlatk
call ZTLLC(ALMD,APHD,TLM0D,DTR,CTPH0,STPH0,TLM,TPH)
wbd0 = wbd -2*dlmd
sbd0 = sbd-dphd
dlon = 2*dlmd
dlat = 2*dphd
csa= sqrt(2.)/2.
p1 = (tlm - wbd0)/dlon
q1 = (tph - sbd0)/dlat
p = p1-int(p1)
q = q1-int(q1)
c iv = 1 + p1
iv = 0 + p1
jv2= 0 + q1
jv = 2*jv2
pp = p*csa+q*csa
qp =-p*csa+q*csa
c print *,p ,q
c print *,pp,qp
s05 = sqrt(2.0) * 0.5
if ( pp.le.s05.and.qp.ge.0) then
ih = iv
jh = jv + 1
else if ( pp.gt.s05.and.qp.ge.0) then
ih = iv
jh = jv + 2
else if ( pp.le.s05.and.qp.lt.0) then
ih = iv
jh = jv
else if ( pp.gt.s05.and.qp.lt.0) then
ih = iv+1
jh = jv+1
endif
if ( ih.le.0.or.ih.gt.IM.or.jh.le.0.or.jh.gt.JM)then
ih = 0
jh = 0
endif
return
end
subroutine ll2ij ( ar7,IM,JM,rlatk,rlonk,imn,jmn)
PARAMETER (PI=3.141592654, s2d=1.4142136)
PARAMETER (DTR=PI/180.0,dx=1./120.)
dimension ar7(7)
tlm0d = AR7(2)
wbd = AR7(3)
dlmd = AR7(4)
tph0d = AR7(5)
sbd = AR7(6)
dphd = AR7(7)
TPH0=TPH0D*DTR
CTPH0=COS(TPH0)
STPH0=SIN(TPH0)
almd = rlonk
aphd = rlatk
call ZTLLC(ALMD,APHD,TLM0D,DTR,CTPH0,STPH0,TLM,TPH)
dmin = 1000000
imn = 0
jmn = 0
do j=1,JM
do i=1,IM
tlon0 = wbd + mod(j+1,2)*dlmd + (i-1)*2*dlmd
tlat0 = sbd + (j-1)*dphd
dist = sqrt((tlm-tlon0)**2+(tph-tlat0)**2)
if ( dist.le.dmin) then
c if ( dist.le.dmin.and.dhgt(i,j).eq.0) then
imn = i
jmn = j
dmin = dist
endif
enddo
enddo
if ( dmin.ge.5.1) then
imn = 0
jmn = 0
endif
c print *,rlatk,rlonk,imn,jmn
return
end
subroutine vrf4 ( xsr,erm,erma,errms,r,cef,n,UND,ix,TO,TF)
dimension TO(ix), TF(ix)
n = 0
erm = 0
erma = 0
errms= 0
XSRO = 0
XSR = 0
YSR = 0
do i=1, ix
if ( TO(i).ne.UND) then
NSRO = NSRO + 1
XSRO = XSRO+TO(i)
endif
if ( TO(i).ne.UND.and.TF(i).ne.UND) then
XSR = XSR + TO(i)
YSR = YSR + TF(i)
erm = erm + ( TF(i)-TO(i))
erma = erma + abs( TF(i)-TO(i))
errms= errms+ ( TF(i)-TO(i)) * ( TF(i)-TO(i))
n = n + 1
endif
enddo
XSR = XSR/n
YSR = YSR/n
erm = erm/n
erma = erma/n
fv = errms
errms= sqrt ( errms/n )
XSRO = XSRO/NSRO
ov = 0
do i=1, ix
if ( TO(i).ne.UND.and.TF(i).ne.UND) then
ov = ov+(TO(i)-XSR)**2.
endif
enddo
cef = 1.-fv/ov
ZXX=0.
ZYY=0.
ZXY=0.
do i=1, ix
if ( TO(i).ne.UND.and.TF(i).ne.UND) then
XT=TO(i)-XSR
YT=TF(i)-YSR
ZXX=ZXX+XT**2
ZYY=ZYY+YT**2
ZXY=ZXY+XT*YT
endif
enddo
R=ZXY/SQRT(ZXX*ZYY)
return
end
c==========================================================================
subroutine oscmbr(INON,rn0,zn0,sgn20
., ar7,IM,JM,fz,sm,ff,aa,NO,NOM,ola,olo,oo,ao,fo,fo0)
c
dimension ar7(7)
dimension fz(IM,JM),sm(im,jm),ff(IM,JM), aa(IM,JM)
dimension ola(NOM),olo(NOM),oo(NOM),ao(NOM),fo(NOM),fo0(NOM)
dimension tola(1000),tolo(1000),toz(1000),too(1000)
c---------------------------------------------------------------
REAL, ALLOCATABLE:: wij(:,:,:),wijo(:,:)
REAL, ALLOCATABLE:: AMX(:,:),XMX(:),YMX(:)
c---------------------------------------------------------------
if ( NO.gt.1000) then
print *," no analizys NOBS?1000"
stop
endif
UND = -99
c
TLM0D= ar7 (2)
WBD = ar7 (3)
DLMD = ar7 (4)
TPH0D= ar7 (5)
SBD = ar7 (6)
DPHD = ar7 (7)
IM1= -wbd/dlmd +1.5
JM1=2*(-sbd/dphd)+1.5
c print *," im1 jm1 ",im1,jm1
IMJM = im*jm - jm/2
KNE=IM
KNW=IM-1
NINC=2*IM-1
dtr = .01745329
ctph0 = cos(tph0d*dtr)
stph0 = sin(tph0d*dtr)
c=======================================================
c 0) bilin. interp at OBS point!
c=======================================================
ao =-99
NOBS = NO
do k=1, NOBS
call hbox ( ar7,IM,JM,ola(k),olo(k),i,j)
c================================================================
htlm = wbd+(i-1)*2*dlmd+mod(j+1,2)*dlmd
htph = sbd+(j-1)*1*dphd
hzet = fz(i,j)
tolo(k) = htlm
tola(k) = htph
toz(k) = hzet
ao (k) = ff(i,j)
fo0(k) = ff(i,j)
c to(NOBS) = oo(k)
c tf(NOBS) = ff(i,j)
c================================================================
c call bilinlist ( und,ar7,IM,JM,fz,ola(k),olo(k),tlm,tph,bzet)
c call bilinlist ( und,ar7,IM,JM,ff,ola(k),olo(k),tlm,tph,bf )
c tolo(k) = tlm
c tola(k) = tph
c toz(k) = bzet
c ao (k) = bf
c fo0(k) = ff(i,j)
71 continue
c tolo(NOBS) = tlm
c tola(NOBS) = tph
enddo
c
500 continue
c stop
allocate (wij(IM,JM,NOBS))
allocate (wijo(NOBS,NOBS))
allocate (AMX(NOBS,NOBS))
allocate (XMX(NOBS))
allocate (YMX(NOBS))
c call vrf3 ( erm,erma,errms,r,n,NOM,UND,1000,TO,TF)
31 format ( " Ao ,n,erm,erma,errms,r",i3,4f8.3)
32 format ( " bilinFG ,n,erm,erma,errms,r",i3,4f8.3)
33 format ( " bilinAN ,n,erm,erma,errms,r",i3,4f8.3)
c) analisys at grid points = backgrmund !
aa = ff
dtr = .01745329
c=======================================================
sgn2= sgn20
do 100 itr=1, 8 !
c=======================================================
rn = rn0-(itr-1)*rn0/10.
zn = zn0-(itr-1)*zn0/10.
c=======================================================
print *,itr,rn0,rn,zn0,zn
c if ( itr.gt.5) then
c rn = 0.5*rn0
c zn = 0.5*zn0
c endif
RN2 = RN*RN
ZN2 = ZN*ZN
c=======================================================
c 1.1) weightingS!
c forecast to observation !
c=======================================================
c half = 0.5 ! 1 ! 0.5
half = 1 ! 0.5
do j=1,JM
tlat = sbd + (j-1)*dphd
do i=1,IM
tlon = wbd + (i-1)*2*dlmd + dlmd * mod(j+1,2)
do k=1,NOBS
tlm = tolo(k)
tph = tola(k)
tz = toz(k)
rlat = 0.5*(tph+tlat)
rx = cos(rlat*dtr)*111*(tlon-tlm)
ry = 111*(tlat-tph)
r2 = rx*rx+ry*ry
z2 = (fz(i,j)-tz)**2
if ( zn2.gt.0) then
zzz=exp(half*-z2/ZN2)
else
zzz = 1
endif
rrr=exp(half*-r2/RN2)
citt..........................
zitt=1
citt..........................
wij (i,j,k)=zzz * rrr * zitt
c if(i.eq.im/2) print *,i,j,k,wij (i,j,k)
enddo
enddo
enddo
do j=1,NOBS
tlmj = tolo(j)
tphj = tola(j)
tzj = toz(j)
do k=1,NOBS
tlmk = tolo(k)
tphk = tola(k)
tzk = toz(k)
rlat = 0.5*(tphk+tphj)
rx = cos(rlat*dtr)*111*(tlmk-tlmj)
ry = 111*(tphk-tphj)
r2 = rx*rx+ry*ry
z2 = (tzj-tzk)**2
if ( zn2.gt.0) then
zzz=exp(half*-z2/ZN2)
else
zzz = 1
endif
rrr=exp(half*-r2/RN2)
czitt................................
zitt=1
czitt................................
wijo(j,k) = zzz * rrr *zitt
enddo
enddo
c interpolation
c analisys of observation (ao) to grid points (aa):
alph = 1 ! 0.65
AMX = 0
nsm = 0
do j=1,JM
do i=1,IM
c if ( sm(i,j).eq.1) goto 600
nsm = nsm +1
rim = sgn2
do k=1,NOBS
rim=rim+wij(i,j,k)
enddo
rim = rim*alph
do k=1,NOBS
amx(k,k) = wij(i,j,k)/rim
XMX(k) = (oo(k)-ao(k))
enddo
c====================================
YMX = matmul(AMX,XMX)
c====================================
do k=1,NOBS
aa(i,j) = aa(i,j)+YMX(k)
enddo
600 continue
enddo
enddo
201 continue
c print *," nsm=",nsm,im*jm
c stop
c====================================
c anlysed values at OBS point!
eps2 = sgn2
c==============
fo = ao ! pazi
c==============
do j=1, NOBS
rim = sgn2
do k=1,NOBS
rim=rim+wijo(j,k)
enddo
rim = rim * alph
XMX(j) = (oo(j)-ao(j))
do k=1,NOBS
d = 0
if( k.eq.j) d = 1
AMX(j,k) = (wijo(j,k)+d*eps2)/rim
enddo
enddo ! of j
c-----------------------------
YMX = matmul(AMX,XMX)
c print *,itr," endof matmull rn zn",rn,zn,nobs
c-----------------------------
do j=1,NOBS
fo(j) = ao(j)+YMX(j)
enddo
c===================
ao = fo ! pazi
c=======================================================
100 continue ! iteration
c=======================================================
c===================
if ( INON.eq.1) then
snp =0
snn = 0
do j=1,jm
do i=1,im
if ( aa(i,j).lt.0) then
snn = snn + aa(i,j)
else
snp = snp + aa(i,j)
endif
enddo
enddo
do j=1,jm
do i=1,im
aa(i,j)=(1+snn/snp)*aa(i,j)
aa(i,j)=max(0.,aa(i,j))
enddo
enddo
if(snp.gt.0) then
c print *," nonnegative",snn/snp
endif
endif
c goto 100
c1) error of analisys value:
c=======================================================
c 100 continue ! iteration
c=======================================================
deallocate (wij,stat=ier)
deallocate (wijo,stat=ier)
deallocate (AMX,stat=ier)
deallocate (XMX,stat=ier)
deallocate (YMX,stat=ier)
return
end
subroutine tocheckll1 (rlon,rlat,GWBD,GSBD,GEBD,GNBD,IRET)
iret = 0
if( rlon.lt.GWBD.or.rlon.gt.GEBD) then
print *," westeastHH stop"
print *,rlat,rlon,GWBD,GEBD
iret = 1
stop
endif
if( rlat.lt.GSBD.or.rlat.gt.GNBD) then
print *," southnorthHH stop"
print *,rlat,rlon,GSBD,GNBD
iret = 1
stop
endif
return
end
subroutine fdiskm ( rla,rlo,dstmn,i1,NN,mst,cla,clo)
dimension cla(mst),clo(mst)
dstmn = 1.e+9
do i=1, NN
defi = rla-cla(i)
dela = rlo-clo(i)
pk = (sind(0.5*defi))**2.
+ +cosd(rla)*cosd(cla(i))*((sind(0.5*dela))**2.)
dst = 6371*2.*asin(sqrt(pk))
if(dst.le.dstmn) then
dstmn = dst
i1 = i
endif
enddo
return
end