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 if (ar7(1).eq.3.) then call b2ll (omn,omx,ar7,iuv,jfile,i9,j9,vl1(1,1,L),vl2(1,1,L)) 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 b2ll (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) cbmreza im=-2.*wbd/dlmd+1.5 jm=-2.*sbd/dphd+1.5 iglobal = 0 if(abs(wbd).eq.180) iglobal = 1 if(iglobal.eq.1) then im = im+2 jm = jm+2 endif wbd1=wbd sbd1=sbd if( iglobal.eq.1) then wbd1=wbd -dlmd !/2. sbd1=sbd -dphd ! *0.5 endif if(iuv.eq.2) then wbd1=wbd+dlmd/2. sbd1=sbd+dphd*0.5 endif 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 if(iglobal.eq.0) then call ZTLLC(X,Y,TLM0D,DTR,CTPH0,STPH0,TLM,TPH) else tlm = X tph = Y endif 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 call Zgcoef (ar7,tlm,tph, kh,ph,qh, kv,pv,qv) ! i0=(tlm-wbd1)/dlmd+1.5 ! j0=(tph-sbd1)/dphd+1.5 p=(tlm-wbd1)/dlmd ph=p-int(p) q=(tph-sbd1)/dphd qh=q-int(q) i0 = int(p)+1 j0 = int(q)+1 kh1= (j0-1)*im+i0 kh2 = kh1 + 1 kh3= kh2+im kh4 = kh3 - 1 ! print *,wbd,wbd1,sbd,sbd1 ! print *,ph,qh,i0,j0 ! print *,tlm,tph ! print *,wbd1+(i0-1)*dlmd,wbd1+(i0-0)*dlmd ! print *,sbd1+(j0-1)*dphd,sbd1+(j0-0)*dphd if(ph.lt.0.or.qh.lt.0) then print *,"negative ph qh",ph,qh stop endif ! stop 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)) if( iglobal.eq.0) then call ZRLTLW (X, Y, v1, v2, tlm0d,dtr,ctph0,stph0 + , vall1(ii), vall2(ii)) else vall1(ii) = v1 vall2(ii) = v2 endif 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 print *,'iglobal',iglobal,' im jm',im,jm c print *,val1(1:im*jm) return end 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