Subversion Repositories lagranto.ecmwf

Rev

Rev 3 | Rev 13 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 3 Rev 9
Line 100... Line 100...
100 
c     Auxiliary variables
100 
c     Auxiliary variables
101 
      real                                   delta,rd
 101 
      real                                   delta,rd
102 
      integer	                             itm,iloop,i,j,k,filo,lalo
 102 
      integer	                             itm,iloop,i,j,k,filo,lalo
103 
      integer                                ierr,stat
 103 
      integer                                ierr,stat
104 
      integer                                cdfid,fid
 104 
      integer                                cdfid,fid
105 
      real	                             tstart,time0,time1,time
 105 
      real	                                 tstart,time0,time1,time
106 
      real                                   reltpos0,reltpos1
 106 
      real                                   reltpos0,reltpos1
107 
      real                                   xind,yind,pind,pp,sp,stagz
 107 
      real                                   xind,yind,pind,pp,sp,stagz
108 
      character*80                           filename,varname
 108 
      character*80                           filename,varname
109 
      integer                                reftmp(6)
 109 
      integer                                reftmp(6)
110 
      character                              ch
 110 
      character                              ch
111 
      real                                   frac,tload
 111 
      real                                   frac,tload
112 
      integer                                itim
 112 
      integer                                itim
113 
      integer                                wstep
 113 
      integer                                wstep
- 
 
 114 
      real                                   x1,y1
- 
 
 115 
 
- 
 
 116 
c     Externals
- 
 
 117 
      real                                   int_index4
- 
 
 118 
      external                               int_index4
114 
 
 119 
 
115 
c     --------------------------------------------------------------------
 120 
c     --------------------------------------------------------------------
116 
c     Start of program, Read parameters
 121 
c     Start of program, Read parameters
117 
c     --------------------------------------------------------------------
 122 
c     --------------------------------------------------------------------
118 
 
 123 
 
Line 253... Line 258...
253 
 
 258 
 
254 
C     Check if the number of levels is too large
 259 
C     Check if the number of levels is too large
255 
      if (nz.gt.nlevmax) goto 993
 260 
      if (nz.gt.nlevmax) goto 993
256 
 
 261 
 
257 
C     Set logical flag for periodic data set (hemispheric or not)
 262 
C     Set logical flag for periodic data set (hemispheric or not)
- 
 
 263 
      hem = 0
258 
      if (per.eq.0.) then
 264 
      if (per.eq.0.) then
259 
         delta=xmax-xmin-360.
 265 
         delta=xmax-xmin-360.
260 
         if (abs(delta+dx).lt.eps) then               ! Program aborts: arrays must be closed
 266 
         if (abs(delta+dx).lt.eps) then               ! Program aborts: arrays must be closed
261 
            goto 992
 267 
            goto 992
262 
        else if (abs(delta).lt.eps) then              ! Periodic and hemispheric
 268 
        else if (abs(delta).lt.eps) then              ! Periodic and hemispheric
Line 372... Line 378...
372 
         traout(i,itim,1) = 0.
 378 
         traout(i,itim,1) = 0.
373 
         traout(i,itim,2) = xx0(i)
 379 
         traout(i,itim,2) = xx0(i)
374 
         traout(i,itim,3) = yy0(i)
 380 
         traout(i,itim,3) = yy0(i)
375 
         traout(i,itim,4) = pp0(i)
 381 
         traout(i,itim,4) = pp0(i)
376 
      enddo
 382 
      enddo
377
383 
 
378 
c     Init the flag and the counter for trajectories leaving the domain
 384 
c     Init the flag and the counter for trajectories leaving the domain
379 
      leftcount=0
 385 
      leftcount=0
380 
      do i=1,ntra
 386 
      do i=1,ntra
381 
         leftflag(i)=0
 387 
         leftflag(i)=0
382 
      enddo
 388 
      enddo
Line 385... Line 391...
385 
      call hhmm2frac(tst,frac)
 391 
      call hhmm2frac(tst,frac)
386 
      tst = frac
 392 
      tst = frac
387 
      call hhmm2frac(ten,frac)
 393 
      call hhmm2frac(ten,frac)
388 
      ten = frac
 394 
      ten = frac
389 
 
 395 
 
- 
 
 396 
c     Check that all starting positions are above topography
- 
 
 397 
      varname = 'P'
- 
 
 398 
      filename = prefix//dat(1)
- 
 
 399 
      call input_open (fid,filename)
- 
 
 400 
      call input_grid
- 
 
 401 
     >    (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
- 
 
 402 
     >     tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
- 
 
 403 
      call input_close(fid)
- 
 
 404 
      do i=1,ntra
- 
 
 405 
 
- 
 
 406 
C       Interpolate surface pressure to actual position (from first input file)
- 
 
 407 
        x1 = xx0(i)
- 
 
 408 
        y1 = yy0(i)
- 
 
 409 
        call get_index4 (xind,yind,pind,x1,y1,1050.,0.,
- 
 
 410 
     >                   p3t1,p3t1,spt1,spt1,3,
- 
 
 411 
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 412 
        sp = int_index4 (spt1,spt1,nx,ny,1,xind,yind,1.,0.,mdv)
- 
 
 413 
 
- 
 
 414 
c       Decide whether to keep the trajectory
- 
 
 415 
        if ( pp0(i).gt.sp ) then
- 
 
 416 
            write(*,'(a30,3f10.2)')
- 
 
 417 
     >               'WARNING: starting point below topography ',
- 
 
 418 
     >               xx0(i),yy0(i),pp0(i)
- 
 
 419 
            leftflag(i) = 1
- 
 
 420 
        endif
- 
 
 421 
 
- 
 
 422 
      enddo
- 
 
 423 
 
- 
 
 424 
 
390 
c     -----------------------------------------------------------------------
 425 
c     -----------------------------------------------------------------------
391 
c     Loop to calculate trajectories
 426 
c     Loop to calculate trajectories
392 
c     -----------------------------------------------------------------------
 427 
c     -----------------------------------------------------------------------
393 
 
 428 
 
394 
c     Write some status information
 429 
c     Write some status information
- 
 
 430 
      print*
395 
      print*,'---- TRAJECTORIES ----------- ---------------------------'
 431 
      print*,'---- TRAJECTORIES ----------- ---------------------------'
396 
      print*
432 
      print*
397 
 
 433 
 
398 
C     Set the time for the first data file (depending on forward/backward mode)
 434 
C     Set the time for the first data file (depending on forward/backward mode)
399 
      if (fbflag.eq.1) then
 435 
      if (fbflag.eq.1) then
Line 428... Line 464...
428 
     >     xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
464 
     >     xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
429 
      call input_grid                                  ! GRID
 465 
      call input_grid                                  ! GRID
430 
     >    (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
 466 
     >    (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
431 
     >     tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
 467 
     >     tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
432 
      call input_close(fid)
 468 
      call input_close(fid)
433 
 
 469
434 
c     Loop over all input files (time step is <timeinc>)
 470 
c     Loop over all input files (time step is <timeinc>)
435 
      do itm=1,numdat-1
 471 
      do itm=1,numdat-1
436 
 
 472 
 
437 
c       Calculate actual and next time
 473 
c       Calculate actual and next time
438 
        time0 = tstart+real(itm-1)*timeinc*fbflag
 474 
        time0 = tstart+real(itm-1)*timeinc*fbflag
Line 494... Line 530...
494 
 
 530 
 
495 
C         Calculate relative time position in the interval timeinc (0=beginning, 1=end)
 531 
C         Calculate relative time position in the interval timeinc (0=beginning, 1=end)
496 
          reltpos0 = ((real(iloop)-1.)*ts)/timeinc
 532 
          reltpos0 = ((real(iloop)-1.)*ts)/timeinc
497 
          reltpos1 = real(iloop)*ts/timeinc
 533 
          reltpos1 = real(iloop)*ts/timeinc
498 
 
 534 
 
499 
C         Initialize counter for domain leaving trajectories
 - 
 
500 
          leftcount=0
 - 
 
501 
 
 - 
 
502 
c         Timestep for all trajectories
 535 
c         Timestep for all trajectories
503 
          do i=1,ntra
 536 
          do i=1,ntra
504 
 
 537 
 
505 
C           Check if trajectory has already left the data domain
 538 
C           Check if trajectory has already left the data domain
506 
            if (leftflag(i).ne.1) then
539 
            if (leftflag(i).ne.1) then
Line 554... Line 587...
554 
          if (abs(delta).lt.eps) then
 587 
          if (abs(delta).lt.eps) then
555 
c          wstep = wstep + abs(ts)
 588 
c          wstep = wstep + abs(ts)
556 
c          if ( mod(wstep,deltout).eq.0 ) then
 589 
c          if ( mod(wstep,deltout).eq.0 ) then
557 
            time = time0+reltpos1*timeinc*fbflag
 590 
            time = time0+reltpos1*timeinc*fbflag
558 
            itim = itim + 1
 591 
            itim = itim + 1
- 
 
 592 
            if ( itim.le.ntim ) then
559 
            do i=1,ntra
 593 
              do i=1,ntra
560 
               call frac2hhmm(time,tload)
 594 
                 call frac2hhmm(time,tload)
561 
               traout(i,itim,1) = tload
 595 
                 traout(i,itim,1) = tload
562 
               traout(i,itim,2) = xx0(i)
 596 
                 traout(i,itim,2) = xx0(i)
563 
               traout(i,itim,3) = yy0(i)
597 
                 traout(i,itim,3) = yy0(i)
564 
               traout(i,itim,4) = pp0(i)
 598 
                 traout(i,itim,4) = pp0(i)
565 
            enddo
 599 
              enddo
- 
 
 600 
            endif
566 
          endif
 601 
          endif
567 
 
 602 
 
568 
        enddo
 603 
        enddo
569 
 
 604 
 
570 
      enddo
 605 
      enddo
571 
 
 606 
 
- 
 
 607 
      print*,(traout(1,1,i),i=1,4)
- 
 
 608 
 
572 
c     Write trajectory file
 609 
c     Write trajectory file
573 
      vars(1)  ='time'
 610 
      vars(1)  ='time'
574 
      vars(2)  ='lon'
 611 
      vars(2)  ='lon'
575 
      vars(3)  ='lat'
 612 
      vars(3)  ='lat'
576 
      vars(4)  ='p'
 613 
      vars(4)  ='p'