Subversion Repositories lagranto.wrf

Rev

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

Rev 19 Rev 21
Line 69... Line 69...
69 
      integer                                jflag           ! Jump flag (if =1 ground-touching trajectories reenter atmosphere)
 69 
      integer                                jflag           ! Jump flag (if =1 ground-touching trajectories reenter atmosphere)
70 
      real                                   wfactor         ! Factor for vertical velocity field
 70 
      real                                   wfactor         ! Factor for vertical velocity field
71 
      character*80                           strname         ! File with start positions
 71 
      character*80                           strname         ! File with start positions
72 
      character*80                           timecheck       ! Either 'yes' or 'no'
 72 
      character*80                           timecheck       ! Either 'yes' or 'no'
73 
      character*10		              	     gridflag        ! "ideal" or "real"
 73 
      character*10		              	     gridflag        ! "ideal" or "real"
- 
 
 74 
      real                                   aglconst        ! constant level above ground
- 
 
 75 
      real                                   isen            ! isentropic level
- 
 
 76 
      character*80                           vertmode        ! Vertical mode (kinematic,aglconst)
74 
 
 77 
 
75 
c     Trajectories
 78 
c     Trajectories
76 
      integer                                ncol            ! Number of columns for insput trajectories
 79 
      integer                                ncol            ! Number of columns for insput trajectories
77 
      real,allocatable, dimension (:,:,:) :: trainp          ! Input start coordinates (ntra,1,ncol)
 80 
      real,allocatable, dimension (:,:,:) :: trainp          ! Input start coordinates (ntra,1,ncol)
78 
      real,allocatable, dimension (:,:,:) :: traout          ! Output trajectories (ntra,ntim,4)
 81 
      real,allocatable, dimension (:,:,:) :: traout          ! Output trajectories (ntra,ntim,4)
Line 88... Line 91...
88 
      real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
 91 
      real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
89 
      real,allocatable, dimension (:)     :: uut0,uut1       ! Zonal wind
 92 
      real,allocatable, dimension (:)     :: uut0,uut1       ! Zonal wind
90 
      real,allocatable, dimension (:)     :: vvt0,vvt1       ! Meridional wind
 93 
      real,allocatable, dimension (:)     :: vvt0,vvt1       ! Meridional wind
91 
      real,allocatable, dimension (:)     :: wwt0,wwt1       ! Vertical wind
 94 
      real,allocatable, dimension (:)     :: wwt0,wwt1       ! Vertical wind
92 
      real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure
95 
      real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure
- 
 
 96 
      real,allocatable, dimension (:)     :: tht0,tht1       ! potential temperature
- 
 
 97 
      real,allocatable, dimension (:)     :: sth0,sth1       ! surface potential temperature
93 
 
 98 
 
94 
c     Grid description
 99 
c     Grid description
95 
      real                                   pollon,pollat   ! Longitude/latitude of pole
 100 
      real                                   pollon,pollat   ! Longitude/latitude of pole
96 
      real                                   ak(nlevmax)     ! Vertical layers and levels
 101 
      real                                   ak(nlevmax)     ! Vertical layers and levels
97 
      real                                   bk(nlevmax)
102 
      real                                   bk(nlevmax)
Line 119... Line 124...
119 
      integer                                wstep
 124 
      integer                                wstep
120 
      character*80                           radunit
 125 
      character*80                           radunit
121 
      real                                   lon1,lat1,lon2,lat2
 126 
      real                                   lon1,lat1,lon2,lat2
122 
      real                                   scale_max,scale_min
 127 
      real                                   scale_max,scale_min
123 
 
 128 
 
- 
 
 129 
c     Externals
- 
 
 130 
      real         int_index4
- 
 
 131 
      external     int_index4
- 
 
 132 
 
124 
c     --------------------------------------------------------------------
 133 
c     --------------------------------------------------------------------
125 
c     Start of program, Read parameters
 134 
c     Start of program, Read parameters
126 
c     --------------------------------------------------------------------
 135 
c     --------------------------------------------------------------------
127 
 
 136 
 
128 
c     Write start message
 137 
c     Write start message
Line 195... Line 204...
195 
      read(9,*) reftime(5)                  ! min
 204 
      read(9,*) reftime(5)                  ! min
196 
      read(9,*) reftime(6)                  ! time range (in min)
 205 
      read(9,*) reftime(6)                  ! time range (in min)
197 
 
 206 
 
198 
c     Read flag for 'no time check'
 207 
c     Read flag for 'no time check'
199 
      read(9,*) timecheck
 208 
      read(9,*) timecheck
- 
 
 209
- 
 
 210 
c     Read flag for trajectories at constant level above ground
- 
 
 211 
      read(9,*) aglconst
- 
 
 212
- 
 
 213 
c     Read flag for trajectories at constant isentropic level
- 
 
 214 
      read(9,*) isen
200 
 
 215 
 
201 
c     Close the input file
 216 
c     Close the input file
202 
      close(9)
 217 
      close(9)
203 
 
 218 
 
- 
 
 219 
c     Set the vertical mode of the trajectory calculation
- 
 
 220 
      vertmode = 'kinematic'
- 
 
 221 
      if ( aglconst.gt.0. ) then
- 
 
 222 
         vertmode = 'aglconst'
- 
 
 223 
      endif
- 
 
 224 
      if ( isen.gt.0. ) then
- 
 
 225 
         vertmode = 'isen'
- 
 
 226 
      endif
- 
 
 227 
 
204 
c     Calculate the number of output time steps
 228 
c     Calculate the number of output time steps
205 
      ntim = abs(reftime(6)/deltout) + 1
 229 
      ntim = abs(reftime(6)/deltout) + 1
206 
 
 230 
 
207 
c     Set the formats of the input and output files
 231 
c     Set the formats of the input and output files
208 
      call mode_tra(inpmode,strname)
 232 
      call mode_tra(inpmode,strname)
Line 232... Line 256...
232 
         print*,'  Input format           : ',inpmode
 256 
         print*,'  Input format           : ',inpmode
233 
      endif
 257 
      endif
234 
      print*,'  Output format          : ',outmode
 258 
      print*,'  Output format          : ',outmode
235 
      print*,'  Periodicity            : ',per
 259 
      print*,'  Periodicity            : ',per
236 
      print*,'  Time check             : ',trim(timecheck)
 260 
      print*,'  Time check             : ',trim(timecheck)
- 
 
 261 
      if ( vertmode.eq.'aglconst' ) then
- 
 
 262 
         print*,'  Z = const AGL          : ',aglconst
- 
 
 263 
      endif
- 
 
 264 
      if ( vertmode.eq.'isen' ) then
- 
 
 265 
         print*,'  isentropic             : ',isen
- 
 
 266 
      endif
- 
 
 267
237 
      print*
 268 
      print*
238 
 
 269 
 
239 
      print*,'---- FIXED NUMERICAL PARAMETERS -------------------------'
 270 
      print*,'---- FIXED NUMERICAL PARAMETERS -------------------------'
240 
      print*
 271 
      print*
241 
      print*,'  Number of iterations   : ',numit
 272 
      print*,'  Number of iterations   : ',numit
Line 307... Line 338...
307 
      allocate(pp0(ntra),stat=stat)
 338 
      allocate(pp0(ntra),stat=stat)
308 
      if (stat.ne.0) print*,'*** error allocating array pp0      ***' ! Pressure
 339 
      if (stat.ne.0) print*,'*** error allocating array pp0      ***' ! Pressure
309 
      allocate(leftflag(ntra),stat=stat)
 340 
      allocate(leftflag(ntra),stat=stat)
310 
      if (stat.ne.0) print*,'*** error allocating array leftflag ***' ! Leaving-domain flag
 341 
      if (stat.ne.0) print*,'*** error allocating array leftflag ***' ! Leaving-domain flag
311 
 
 342 
 
- 
 
 343 
c     Get memory for potential temperature (isentropic trajectories)
- 
 
 344 
      if ( vertmode.eq.'isen' ) then
- 
 
 345 
          allocate(tht0(nx*ny*nz),stat=stat)
- 
 
 346 
          if (stat.ne.0) print*,'*** error allocating array tht0 ***'
- 
 
 347 
          allocate(tht1(nx*ny*nz),stat=stat)
- 
 
 348 
          if (stat.ne.0) print*,'*** error allocating array tht1 ***'
- 
 
 349 
          allocate(sth0(nx*ny),stat=stat)
- 
 
 350 
          if (stat.ne.0) print*,'*** error allocating array sth0 ***'
- 
 
 351 
          allocate(sth1(nx*ny),stat=stat)
- 
 
 352 
          if (stat.ne.0) print*,'*** error allocating array sth1 ***'
- 
 
 353 
      endif
- 
 
 354 
 
312 
c     Write some status information
 355 
c     Write some status information
313 
      print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
 356 
      print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
314 
      print*
 357 
      print*
315 
      print*,'  xmin,xmax     : ',xmin,xmax
 358 
      print*,'  xmin,xmax     : ',xmin,xmax
316 
      print*,'  ymin,ymax     : ',ymin,ymax
 359 
      print*,'  ymin,ymax     : ',ymin,ymax
Line 437... Line 480...
437 
     >        ( reftime(5).ne.reftmp(5) ) )
 480 
     >        ( reftime(5).ne.reftmp(5) ) )
438 
     >   then
 481 
     >   then
439 
            print*,' WARNING: Inconsistent reference times'
 482 
            print*,' WARNING: Inconsistent reference times'
440 
            write(*,'(5i8)') (reftime(i),i=1,5)
 483 
            write(*,'(5i8)') (reftime(i),i=1,5)
441 
            write(*,'(5i8)') (reftmp (i),i=1,5)
 484 
            write(*,'(5i8)') (reftmp (i),i=1,5)
442 
            print*,'Enter a key to proceed...'
 485 
            print*,' Enter a key to proceed...'
443 
            stop
 486 
            stop
444 
         endif
 487 
         endif
445 
      endif
 488 
      endif
446 
 
 489 
 
447 
c     Transform start coordinates from index space to WRF grid
 490 
c     Transform start coordinates from index space to WRF grid
Line 462... Line 505...
462 
      print*,'                (hour)  :',reftime(4)
 505 
      print*,'                (hour)  :',reftime(4)
463 
      print*,'                (min)   :',reftime(5)
 506 
      print*,'                (min)   :',reftime(5)
464 
      print*,' Time range             :',reftime(6),' min'
 507 
      print*,' Time range             :',reftime(6),' min'
465 
      print*
 508 
      print*
466 
 
 509 
 
467 
C     Save starting positions
510 
C     Save starting positions - vertical position will be corrected below for <aglconst> and <isen>
468 
      itim = 1
 511 
      itim = 1
469 
      do i=1,ntra
 512 
      do i=1,ntra
470 
         traout(i,itim,1) = 0.
 513 
         traout(i,itim,1) = 0.
471 
         traout(i,itim,2) = xx0(i)
 514 
         traout(i,itim,2) = xx0(i)
472 
         traout(i,itim,3) = yy0(i)
 515 
         traout(i,itim,3) = yy0(i)
Line 483... Line 526...
483 
      call hhmm2frac(tst,frac)
 526 
      call hhmm2frac(tst,frac)
484 
      tst = frac
 527 
      tst = frac
485 
      call hhmm2frac(ten,frac)
 528 
      call hhmm2frac(ten,frac)
486 
      ten = frac
 529 
      ten = frac
487 
 
 530 
 
- 
 
 531 
c     Adjust the starting heights in mode <aglconst>
- 
 
 532 
      if ( vertmode.eq.'aglconst' ) then
- 
 
 533 
        filename = prefix//dat(1)
- 
 
 534 
        call input_open (fid,filename)
- 
 
 535 
	      varname='geopot'				! geopot.height
- 
 
 536 
	      call input_var_wrf(fid,varname,p3t0,nx,ny,nz)
- 
 
 537 
	      varname='geopots'				! geopot.height at surface
- 
 
 538 
	      call input_var_wrf(fid,varname,spt0,nx,ny,1)
- 
 
 539 
        call input_close(fid)
- 
 
 540 
        filename = prefix//dat(2)
- 
 
 541 
        call input_open (fid,filename)
- 
 
 542 
	      varname='geopot'				! geopot.height
- 
 
 543 
	      call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
- 
 
 544 
	      varname='geopots'				! geopot.height at surface
- 
 
 545 
	      call input_var_wrf(fid,varname,spt1,nx,ny,1)
- 
 
 546 
        call input_close(fid)
- 
 
 547 
        do i=1,ntra
- 
 
 548 
           call get_index4 (xind,yind,pind,xx0(i),yy0(i),0.,reltpos0,
- 
 
 549 
     >                      p3t0,p3t1,spt0,spt1,3,
- 
 
 550 
     >                      nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 551 
           sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos0,mdv)
- 
 
 552 
           traout(i,1,4) = sp + aglconst
- 
 
 553 
        enddo
- 
 
 554 
      endif
- 
 
 555
- 
 
 556 
c     Adjust the starting heights in mode <isen>
- 
 
 557 
      if ( vertmode.eq.'isen' ) then
- 
 
 558 
        filename = prefix//dat(1)
- 
 
 559 
        call input_open (fid,filename)
- 
 
 560 
	      varname='geopot'				! geopot.height
- 
 
 561 
	      call input_var_wrf(fid,varname,p3t0,nx,ny,nz)
- 
 
 562 
	      varname='geopots'				! geopot.height at surface
- 
 
 563 
	      call input_var_wrf(fid,varname,spt0,nx,ny,1)
- 
 
 564 
          varname='TH'					! TH
- 
 
 565 
          call input_var_wrf(fid,varname,tht0,nx,ny,nz)
- 
 
 566 
          do i=1,nx*ny
- 
 
 567 
              sth0(i) = tht0(i)
- 
 
 568 
          enddo
- 
 
 569 
        call input_close(fid)
- 
 
 570 
        filename = prefix//dat(2)
- 
 
 571 
        call input_open (fid,filename)
- 
 
 572 
	      varname='geopot'				! geopot.height
- 
 
 573 
	      call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
- 
 
 574 
	      varname='geopots'				! geopot.height at surface
- 
 
 575 
	      call input_var_wrf(fid,varname,spt1,nx,ny,1)
- 
 
 576 
          varname='TH'					! TH
- 
 
 577 
          call input_var_wrf(fid,varname,tht1,nx,ny,nz)
- 
 
 578 
          do i=1,nx*ny
- 
 
 579 
              sth1(i) = tht1(i)
- 
 
 580 
          enddo
- 
 
 581 
        call input_close(fid)
- 
 
 582 
        do i=1,ntra
- 
 
 583
- 
 
 584 
          call get_index4 (xind,yind,pind,xx0(i),yy0(i),isen,reltpos0,
- 
 
 585 
     >                    tht0,tht1,sth0,sth1,1,
- 
 
 586 
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 587 
          traout(i,1,4) = int_index4
- 
 
 588 
     >          (p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 589 
        enddo
- 
 
 590 
      endif
- 
 
 591 
 
488 
c     -----------------------------------------------------------------------
 592 
c     -----------------------------------------------------------------------
489 
c     Loop to calculate trajectories
 593 
c     Loop to calculate trajectories
490 
c     -----------------------------------------------------------------------
 594 
c     -----------------------------------------------------------------------
491 
 
 595 
 
492 
c     Write some status information
 596 
c     Write some status information
Line 526... Line 630...
526 
	  call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
 630 
	  call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
527 
 
 631 
 
528 
	  varname='geopots'				! geopot.height at surface
 632 
	  varname='geopots'				! geopot.height at surface
529 
	  call input_var_wrf(fid,varname,spt1,nx,ny,1)
 633 
	  call input_var_wrf(fid,varname,spt1,nx,ny,1)
530 
 
 634 
 
- 
 
 635 
      if (vertmode.eq.'isen') then
- 
 
 636 
           varname='TH'					! TH
- 
 
 637 
 	       call input_var_wrf(fid,varname,tht1,nx,ny,nz)
- 
 
 638 
           do i=1,nx*ny
- 
 
 639 
              sth1(i) = tht1(i)
- 
 
 640 
           enddo
- 
 
 641 
      endif
- 
 
 642 
 
531 
      call input_close(fid)
 643 
      call input_close(fid)
532 
 
 644 
 
533 
c     Loop over all input files (time step is <timeinc>)
 645 
c     Loop over all input files (time step is <timeinc>)
534 
      do itm=1,numdat-1
 646 
      do itm=1,numdat-1
535 
 
 647 
 
536 
c       Calculate actual and next time
 648 
c       Calculate actual and next time
537 
        time0 = tstart+real(itm-1)*timeinc*fbflag
 649 
        time0 = tstart+real(itm-1)*timeinc*fbflag
538 
        time1 = time0+timeinc*fbflag
 650 
        time1 = time0+timeinc*fbflag
539 
 
 651 
 
540 
 
 - 
 
541 
c       Copy old velocities and pressure fields to new ones
652 
c       Copy old velocities and pressure fields to new ones
542 
        do i=1,nx*ny*nz
 653 
        do i=1,nx*ny*nz
543 
           uut0(i)=uut1(i)
 654 
           uut0(i)=uut1(i)
544 
           vvt0(i)=vvt1(i)
 655 
           vvt0(i)=vvt1(i)
545 
           wwt0(i)=wwt1(i)
 656 
           wwt0(i)=wwt1(i)
Line 547... Line 658...
547 
        enddo
 658 
        enddo
548 
        do i=1,nx*ny
 659 
        do i=1,nx*ny
549 
           spt0(i)=spt1(i)
 660 
           spt0(i)=spt1(i)
550 
        enddo
 661 
        enddo
551 
 
 662 
 
- 
 
 663 
c       Copy potential temperature
- 
 
 664 
        if ( vertmode.eq.'isen' ) then
- 
 
 665 
           do i=1,nx*ny*nz
- 
 
 666 
              tht0(i)=tht1(i)
- 
 
 667 
           enddo
- 
 
 668 
           do i=1,nx*ny
- 
 
 669 
             sth0(i)=sth1(i)
- 
 
 670 
           enddo
- 
 
 671 
        endif
- 
 
 672 
 
552 
c       Read wind fields and surface pressure at next time
 673 
c       Read wind fields and surface pressure at next time
553 
        filename = prefix//dat(itm+1)
 674 
        filename = prefix//dat(itm+1)
554 
 
 675 
 
555 
        call frac2hhmm(time1,tload)
 676 
        call frac2hhmm(time1,tload)
556 
        write(*,'(a16,a20,f7.2)') '  (file,time) : ',
 677 
        write(*,'(a16,a20,f7.2)') '  (file,time) : ',
Line 570... Line 691...
570 
	    varname='geopot'				! geopot.height
 691 
	    varname='geopot'				! geopot.height
571 
	    call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
 692 
	    call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
572 
 
 693 
 
573 
	    varname='geopots'				! geopot.height
 694 
	    varname='geopots'				! geopot.height
574 
	    call input_var_wrf(fid,varname,spt1,nx,ny,1)
 695 
	    call input_var_wrf(fid,varname,spt1,nx,ny,1)
- 
 
 696
- 
 
 697 
        if (vertmode.eq.'isen') then
- 
 
 698 
           varname='TH'					! TH
- 
 
 699 
 	       call input_var_wrf(fid,varname,tht1,nx,ny,nz)
- 
 
 700 
           do i=1,nx*ny
- 
 
 701 
              sth1(i) = tht1(i)
- 
 
 702 
           enddo
- 
 
 703 
        endif
575 
 
 704
576 
        call input_close(fid)
 705 
        call input_close(fid)
577
706
578 
C       Determine the first and last loop indices
 707 
C       Determine the first and last loop indices
579 
        if (numdat.eq.2) then
 708 
        if (numdat.eq.2) then
580 
          filo = nint(tst/ts)+1
 709 
          filo = nint(tst/ts)+1
Line 599... Line 728...
599 
 
 728 
 
600 
C         Initialize counter for domain leaving trajectories
 729 
C         Initialize counter for domain leaving trajectories
601 
          leftcount=0
 730 
          leftcount=0
602 
 
 731 
 
603 
c         Timestep for all trajectories
 732 
c         Timestep for all trajectories
604 
          do i=1,ntra
 733 
          do i=1,ntra
605 
 
 734 
 
606 
C           Check if trajectory has already left the data domain
 735 
C           Check if trajectory has already left the data domain
607 
            if (leftflag(i).ne.1) then
736 
            if (leftflag(i).ne.1) then
608 
 
 737 
 
609 
c             Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
 738 
c             Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
- 
 
 739 
              if ( vertmode.eq.'kinematic' ) then
610 
              call euler(
 740 
                 call euler_kinematic(
611 
     >               xx1,yy1,pp1,leftflag(i),
 741 
     >               xx1,yy1,pp1,leftflag(i),
612 
     >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
 742 
     >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
613 
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
 743 
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
614 
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
 744 
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
615 
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy)
 745 
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy)
- 
 
 746 
              elseif ( vertmode.eq.'aglconst' ) then
- 
 
 747 
                 call euler_aglconst(
- 
 
 748 
     >               xx1,yy1,pp1,leftflag(i),
- 
 
 749 
     >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
- 
 
 750 
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
- 
 
 751 
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,
- 
 
 752 
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,aglconst)
- 
 
 753 
              elseif ( vertmode.eq.'isen' ) then
- 
 
 754 
                 call euler_isen(
- 
 
 755 
     >               xx1,yy1,pp1,leftflag(i),
- 
 
 756 
     >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
- 
 
 757 
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
- 
 
 758 
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,
- 
 
 759 
     >               tht0,tht1,sth0,sth1,
- 
 
 760 
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,isen)
- 
 
 761
- 
 
 762 
              endif
616 
 
 763 
 
617 
c             Update trajectory position, or increase number of trajectories leaving domain
 764 
c             Update trajectory position, or increase number of trajectories leaving domain
618 
              if (leftflag(i).eq.1) then
 765 
              if (leftflag(i).eq.1) then
619 
                leftcount=leftcount+1
 766 
                leftcount=leftcount+1
620 
                if ( leftcount.lt.10 ) then
 767 
                if ( leftcount.lt.10 ) then
Line 721... Line 868...
721 
C     numit	 input	number of iterations
 868 
C     numit	 input	number of iterations
722 
C     jump	 input  flag (=1 trajectories don't enter the ground)
 869 
C     jump	 input  flag (=1 trajectories don't enter the ground)
723 
C     left	 output	flag (=1 if trajectory leaves data domain)
 870 
C     left	 output	flag (=1 if trajectory leaves data domain)
724 
 
 871 
 
725 
c     -------------------------------------------------------------------
 872 
c     -------------------------------------------------------------------
726 
c     Iterative Euler time step
 873 
c     Iterative Euler time step - kinematic trajectory
727 
c     -------------------------------------------------------------------
 874 
c     -------------------------------------------------------------------
728 
 
 875 
 
- 
 
 876 
      subroutine euler_kinematic
729 
      subroutine euler(x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
 877 
     >          (x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
730 
     >                 deltat,numit,jump,mdv,wfactor,fbflag,
 878 
     >           deltat,numit,jump,mdv,wfactor,fbflag,
731 
     >		           spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
 879 
     >		     spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
732 
     >                 xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy)
 880 
     >           xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy)
733 
 
 881 
 
734 
      implicit none
 882 
      implicit none
735 
 
 883 
 
736 
c     Flag for test mode
 884 
c     Flag for test mode
737 
      integer      test
 885 
      integer      test
Line 863... Line 1011...
863 
 
 1011 
 
864 
      return
 1012 
      return
865 
 
 1013 
 
866 
      end
 1014 
      end
867 
 
 1015 
 
- 
 
 1016 
c     -------------------------------------------------------------------
- 
 
 1017 
c     Iterative Euler time step - constant level above ground
- 
 
 1018 
c     -------------------------------------------------------------------
- 
 
 1019 
 
- 
 
 1020 
      subroutine euler_aglconst
- 
 
 1021 
     >          (x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
- 
 
 1022 
     >           deltat,numit,jump,mdv,wfactor,fbflag,
- 
 
 1023 
     >		     spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
- 
 
 1024 
     >           xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,aglconst)
- 
 
 1025 
 
- 
 
 1026 
      implicit none
- 
 
 1027 
 
- 
 
 1028 
c     Flag for test mode
- 
 
 1029 
      integer      test
- 
 
 1030 
      parameter    (test=0)
- 
 
 1031 
 
- 
 
 1032 
c     Declaration of subroutine parameters
- 
 
 1033 
      integer      nx,ny,nz
- 
 
 1034 
      real         x1,y1,p1
- 
 
 1035 
      integer      left
- 
 
 1036 
      real	       x0,y0,p0
- 
 
 1037 
      real         reltpos0,reltpos1
- 
 
 1038 
      real   	   deltat
- 
 
 1039 
      integer      numit
- 
 
 1040 
      integer      jump
- 
 
 1041 
      real         wfactor
- 
 
 1042 
      integer      fbflag
- 
 
 1043 
      real     	   spt0(nx*ny)   ,spt1(nx*ny)
- 
 
 1044 
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
- 
 
 1045 
      real 	       vvt0(nx*ny*nz),vvt1(nx*ny*nz)
- 
 
 1046 
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
- 
 
 1047 
      real         xmin,ymin,dx,dy
- 
 
 1048 
      integer      per
- 
 
 1049 
      integer      hem
- 
 
 1050 
      real         mdv
- 
 
 1051 
      real         mpx(nx*ny),mpy(nx*ny)
- 
 
 1052 
      real         aglconst
- 
 
 1053 
 
- 
 
 1054 
c     Auxiliary variables
- 
 
 1055 
      real         xmax,ymax
- 
 
 1056 
      real	       xind,yind,pind
- 
 
 1057 
 
- 
 
 1058 
      real	       u0,v0,u1,v1,u,v,w,sp
- 
 
 1059 
      integer	   icount
- 
 
 1060 
      character    ch
- 
 
 1061 
      real         mpsc_x,mpsc_y
- 
 
 1062 
 
- 
 
 1063 
c     Externals
- 
 
 1064 
      real         int_index4
- 
 
 1065 
      external     int_index4
- 
 
 1066
- 
 
 1067 
c     Reset the flag for domain-leaving
- 
 
 1068 
      left=0
- 
 
 1069 
 
- 
 
 1070 
c     Set the esat-north boundary of the domain
- 
 
 1071 
      xmax = xmin+real(nx-1)*dx
- 
 
 1072 
      ymax = ymin+real(ny-1)*dy
- 
 
 1073
- 
 
 1074 
C     Set parcel height to aglconst above topography
- 
 
 1075 
      call get_index4 (xind,yind,pind,x1,y1,0.,reltpos0,
- 
 
 1076 
     >                 p3d0,p3d1,spt0,spt1,3,
- 
 
 1077 
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 1078 
      sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos0,mdv)
- 
 
 1079 
      p0 = sp + aglconst
- 
 
 1080 
 
- 
 
 1081 
C     Interpolate wind fields to starting position (x0,y0,p0)
- 
 
 1082 
      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
- 
 
 1083 
     >                 p3d0,p3d1,spt0,spt1,3,
- 
 
 1084 
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 1085 
      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1086 
      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1087 
 
- 
 
 1088 
C     For first iteration take ending position equal to starting position
- 
 
 1089 
      x1 = x0
- 
 
 1090 
      y1 = y0
- 
 
 1091 
      p1 = p0
- 
 
 1092
- 
 
 1093 
C     Iterative calculation of new position
- 
 
 1094 
      do icount=1,numit
- 
 
 1095 
 
- 
 
 1096 
C        Calculate new winds for advection
- 
 
 1097 
         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
- 
 
 1098 
     >                    p3d0,p3d1,spt0,spt1,3,
- 
 
 1099 
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 1100 
         u1 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
- 
 
 1101 
         v1 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
- 
 
 1102
- 
 
 1103 
c        Get the new velocity in between
- 
 
 1104 
         u=(u0+u1)/2.
- 
 
 1105 
         v=(v0+v1)/2.
- 
 
 1106
- 
 
 1107 
c        Get the mapping scale factors for this position
- 
 
 1108 
         mpsc_x = int_index4 (mpx,mpx,nx,ny,1,xind,yind,1.,reltpos1,mdv)
- 
 
 1109 
         mpsc_y = int_index4 (mpy,mpy,nx,ny,1,xind,yind,1.,reltpos1,mdv)
- 
 
 1110 
 
- 
 
 1111 
C        Calculate new positions (adapted for cartesian grid)
- 
 
 1112 
         x1 = x0 + fbflag * deltat * u * dx/mpsc_x
- 
 
 1113 
         y1 = y0 + fbflag * deltat * v * dy/mpsc_y
- 
 
 1114
- 
 
 1115 
c        Write the update positions for tests
- 
 
 1116 
         if ( (icount.eq.3).and.(test.eq.1) ) then
- 
 
 1117 
            write(*,'(10f10.2)') x0,u,x1-x0,p0,w,p1-p0
- 
 
 1118 
         endif
- 
 
 1119 
 
- 
 
 1120 
C       Set trajectory height to aglconst above topography
- 
 
 1121 
        call get_index4 (xind,yind,pind,x1,y1,0.,reltpos1,
- 
 
 1122 
     >                   p3d0,p3d1,spt0,spt1,3,
- 
 
 1123 
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 1124 
        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
- 
 
 1125 
        p1 = sp + aglconst
- 
 
 1126 
 
- 
 
 1127 
c       Handle trajectories which cross the lower boundary (jump flag)
- 
 
 1128 
        if ((jump.eq.1).and.(p1.lt.sp)) then
- 
 
 1129 
            print*,'jump'
- 
 
 1130 
            p1=sp+10.
- 
 
 1131 
        endif
- 
 
 1132 
 
- 
 
 1133 
c       Handle periodioc boundaries for 'ideal' mode
- 
 
 1134 
        if ( per.eq.1 ) then
- 
 
 1135 
	        if ( x1.gt.xmax ) x1=x1-xmax
- 
 
 1136 
	        if ( x1.lt.xmin ) x1=x1+xmax
- 
 
 1137 
        endif
- 
 
 1138
- 
 
 1139 
C       Check if trajectory leaves data domain
- 
 
 1140 
        if ( (x1.lt.xmin).or.(x1.gt.xmax).or.
- 
 
 1141 
     >       (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.lt.sp) )   ! geopot : .lt. ; pressure: .gt.
- 
 
 1142 
     >  then
- 
 
 1143 
          left=1
- 
 
 1144 
          goto 100
- 
 
 1145 
        endif
- 
 
 1146 
 
- 
 
 1147 
      enddo
- 
 
 1148 
 
- 
 
 1149 
c     Exit point for subroutine
- 
 
 1150 
 100  continue
- 
 
 1151 
 
- 
 
 1152 
      return
- 
 
 1153 
 
- 
 
 1154 
      end
- 
 
 1155
- 
 
 1156 
c     -------------------------------------------------------------------
- 
 
 1157 
c     Iterative Euler time step - isentropic
- 
 
 1158 
c     -------------------------------------------------------------------
- 
 
 1159 
 
- 
 
 1160 
      subroutine euler_isen
- 
 
 1161 
     >          (x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
- 
 
 1162 
     >           deltat,numit,jump,mdv,wfactor,fbflag,
- 
 
 1163 
     >		     spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
- 
 
 1164 
     >           tht0,tht1,sth0,sth1,
- 
 
 1165 
     >           xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,isen)
- 
 
 1166 
 
- 
 
 1167 
      implicit none
- 
 
 1168 
 
- 
 
 1169 
c     Flag for test mode
- 
 
 1170 
      integer      test
- 
 
 1171 
      parameter    (test=0)
- 
 
 1172 
 
- 
 
 1173 
c     Declaration of subroutine parameters
- 
 
 1174 
      integer      nx,ny,nz
- 
 
 1175 
      real         x1,y1,p1
- 
 
 1176 
      integer      left
- 
 
 1177 
      real	       x0,y0,p0
- 
 
 1178 
      real         reltpos0,reltpos1
- 
 
 1179 
      real   	   deltat
- 
 
 1180 
      integer      numit
- 
 
 1181 
      integer      jump
- 
 
 1182 
      real         wfactor
- 
 
 1183 
      integer      fbflag
- 
 
 1184 
      real     	   spt0(nx*ny)   ,spt1(nx*ny)
- 
 
 1185 
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
- 
 
 1186 
      real 	       vvt0(nx*ny*nz),vvt1(nx*ny*nz)
- 
 
 1187 
      real         tht0(nx*ny*nz),tht1(nx*ny*nz),sth0(nx*ny),sth1(nx*ny)
- 
 
 1188 
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
- 
 
 1189 
      real         xmin,ymin,dx,dy
- 
 
 1190 
      integer      per
- 
 
 1191 
      integer      hem
- 
 
 1192 
      real         mdv
- 
 
 1193 
      real         mpx(nx*ny),mpy(nx*ny)
- 
 
 1194 
      real         isen
- 
 
 1195 
 
- 
 
 1196 
c     Auxiliary variables
- 
 
 1197 
      real         xmax,ymax
- 
 
 1198 
      real	       xind,yind,pind
- 
 
 1199 
 
- 
 
 1200 
      real	       u0,v0,u1,v1,u,v,w,sp
- 
 
 1201 
      integer	   icount
- 
 
 1202 
      character    ch
- 
 
 1203 
      real         mpsc_x,mpsc_y
- 
 
 1204 
 
- 
 
 1205 
c     Externals
- 
 
 1206 
      real         int_index4
- 
 
 1207 
      external     int_index4
- 
 
 1208
- 
 
 1209 
c     Reset the flag for domain-leaving
- 
 
 1210 
      left=0
- 
 
 1211 
 
- 
 
 1212 
c     Set the esat-north boundary of the domain
- 
 
 1213 
      xmax = xmin+real(nx-1)*dx
- 
 
 1214 
      ymax = ymin+real(ny-1)*dy
- 
 
 1215 
 
- 
 
 1216 
C     Interpolate wind fields to starting position (x0,y0,p0)
- 
 
 1217 
      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
- 
 
 1218 
     >                 p3d0,p3d1,spt0,spt1,3,
- 
 
 1219 
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 1220 
      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1221 
      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1222 
 
- 
 
 1223 
C     For first iteration take ending position equal to starting position
- 
 
 1224 
      x1 = x0
- 
 
 1225 
      y1 = y0
- 
 
 1226 
      p1 = p0
- 
 
 1227
- 
 
 1228 
C     Iterative calculation of new position
- 
 
 1229 
      do icount=1,numit
- 
 
 1230 
 
- 
 
 1231 
C        Calculate new winds for advection
- 
 
 1232 
         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
- 
 
 1233 
     >                    p3d0,p3d1,spt0,spt1,3,
- 
 
 1234 
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 1235 
         u1 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
- 
 
 1236 
         v1 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
- 
 
 1237
- 
 
 1238 
c        Get the new velocity in between
- 
 
 1239 
         u=(u0+u1)/2.
- 
 
 1240 
         v=(v0+v1)/2.
- 
 
 1241
- 
 
 1242 
c        Get the mapping scale factors for this position
- 
 
 1243 
         mpsc_x = int_index4 (mpx,mpx,nx,ny,1,xind,yind,1.,reltpos1,mdv)
- 
 
 1244 
         mpsc_y = int_index4 (mpy,mpy,nx,ny,1,xind,yind,1.,reltpos1,mdv)
- 
 
 1245 
 
- 
 
 1246 
C        Calculate new positions (adapted for cartesian grid)
- 
 
 1247 
         x1 = x0 + fbflag * deltat * u * dx/mpsc_x
- 
 
 1248 
         y1 = y0 + fbflag * deltat * v * dy/mpsc_y
- 
 
 1249
- 
 
 1250 
c        Write the update positions for tests
- 
 
 1251 
         if ( (icount.eq.3).and.(test.eq.1) ) then
- 
 
 1252 
            write(*,'(10f10.2)') x0,u,x1-x0,p0,w,p1-p0
- 
 
 1253 
         endif
- 
 
 1254 
 
- 
 
 1255 
c        Set trajectory height to isentropic level
- 
 
 1256 
         call get_index4 (xind,yind,pind,x1,y1,isen,reltpos1,
- 
 
 1257 
     >                    tht0,tht1,sth0,sth1,1,
- 
 
 1258 
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 1259 
         p1 = int_index4
- 
 
 1260 
     >          (p3d0,p3d1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
- 
 
 1261 
 
- 
 
 1262 
C       Interpolate surface geop.height to actual position
- 
 
 1263 
        call get_index4 (xind,yind,pind,x1,y1,0.,reltpos1,
- 
 
 1264 
     >                   p3d0,p3d1,spt0,spt1,3,
- 
 
 1265 
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
- 
 
 1266 
        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
- 
 
 1267 
 
- 
 
 1268 
c       Handle trajectories which cross the lower boundary (jump flag)
- 
 
 1269 
        if ((jump.eq.1).and.(p1.lt.sp)) then
- 
 
 1270 
            print*,'jump'
- 
 
 1271 
            p1=sp+10.
- 
 
 1272 
        endif
- 
 
 1273 
 
- 
 
 1274 
c       Handle periodioc boundaries for 'ideal' mode
- 
 
 1275 
        if ( per.eq.1 ) then
- 
 
 1276 
	        if ( x1.gt.xmax ) x1=x1-xmax
- 
 
 1277 
	        if ( x1.lt.xmin ) x1=x1+xmax
- 
 
 1278 
        endif
- 
 
 1279
- 
 
 1280 
C       Check if trajectory leaves data domain
- 
 
 1281 
        if ( (x1.lt.xmin).or.(x1.gt.xmax).or.
- 
 
 1282 
     >       (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.lt.sp) )   ! geopot : .lt. ; pressure: .gt.
- 
 
 1283 
     >  then
- 
 
 1284 
          left=1
- 
 
 1285 
          goto 100
- 
 
 1286 
        endif
- 
 
 1287 
 
- 
 
 1288 
      enddo
- 
 
 1289 
 
- 
 
 1290 
c     Exit point for subroutine
- 
 
 1291 
 100  continue
- 
 
 1292 
 
- 
 
 1293 
      return
- 
 
 1294 
 
- 
 
 1295 
      end
868 
 
 1296 
 
869 
c     ----------------------------------------------------------------------
 1297 
c     ----------------------------------------------------------------------
870 
c     Get spherical distance between lat/lon points
 1298 
c     Get spherical distance between lat/lon points
871 
c     ----------------------------------------------------------------------
 1299 
c     ----------------------------------------------------------------------
872 
 
 1300