WebSVN – lagranto.ecmwf – Diff – /trunk/caltra/caltra.f

 c     Input and output format for trajectories (see iotra.f)
       integer   inpmode
       integer   outmode
-c     Filename prefix (typically 'P')
+c     Filename prefix for primary and secondary files (typically 'P' and 'S')
       character*1 prefix
       parameter   (prefix='P')
+      character*1 srefix
+      parameter   (srefix='S')
+c     Physical constants - needed to compute potential temperature
+      real      rdcp,tzero
+      data      rdcp,tzero /0.286,273.15/
 c     --------------------------------------------------------------------
 c     Declaration of variables
 c     --------------------------------------------------------------------
       integer                                deltout         ! Output time interval (in minutes)
       integer                                jflag           ! Jump flag (if =1 ground-touching trajectories reenter atmosphere)
       real                                   wfactor         ! Factor for vertical velocity field
       character*80                           strname         ! File with start positions
       character*80                           timecheck       ! Either 'yes' or 'no'
+      character*80                           isen            ! Isentropic trajectories ('yes' or 'no')
+      integer                                thons           ! Isentropic mode: is TH availanle on S
+      character*80                           modlev          ! 2D (model level) trajectories ('yes' or 'no')
 c     Trajectories
       integer                                ncol            ! Number of columns for insput trajectories
       real,allocatable, dimension (:,:,:) :: trainp          ! Input start coordinates (ntra,1,ncol)
       real,allocatable, dimension (:,:,:) :: traout          ! Output trajectories (ntra,ntim,4)
       real,allocatable, dimension (:)     :: xx0,yy0,pp0     ! Position of air parcels
       integer,allocatable, dimension (:)  :: leftflag        ! Flag for domain-leaving
       real                                   xx1,yy1,pp1     ! Updated position of air parcel
       integer                                leftcount       ! Number of domain leaving trajectories
       integer                                ntim            ! Number of output time steps
+      real,allocatable, dimension (:)     :: theta           ! Potential temperature for isentropic trajectories
+      real,allocatable, dimension (:)     :: zindex          ! Vertical index for model-level (2D) trajectories
 c     Meteorological fields
       real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
       real,allocatable, dimension (:)     :: uut0,uut1       ! Zonal wind
       real,allocatable, dimension (:)     :: vvt0,vvt1       ! Meridional wind
       real,allocatable, dimension (:)     :: wwt0,wwt1       ! Vertical wind
       real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure
+      real,allocatable, dimension (:)     :: tht0,tht1       ! 3d potential temperature
+      real,allocatable, dimension (:)     :: sth0,sth1       ! Surface potential temperature
 c     Grid description
       real                                   pollon,pollat   ! Longitude/latitude of pole
       real                                   ak(nlevmax)     ! Vertical layers and levels
       real                                   bk(nlevmax)
       integer                                reftmp(6)
       character                              ch
       real                                   frac,tload
       integer                                itim
       integer                                wstep
-      real                                   x1,y1
+      real                                   x1,y1,p1
+      real                                   thetamin,thetamax
+      real                                   zindexmin,zindexmax
 c     Externals
       real                                   int_index4
       external                               int_index4
       read(9,*) reftime(6)                  ! time range (in min)
 c     Read flag for 'no time check'
       read(9,*) timecheck
+c     Read flag for isentropic trajectories
+      read(9,*) isen, thons
+c     Read flag for model-level trajectories (2D mode)
+      read(9,*) modlev
 c     Close the input file
       close(9)
 c     Calculate the number of output time steps
       ntim = abs(reftime(6)/deltout) + 1
          print*,'  Input format           : ',inpmode
       endif
       print*,'  Output format          : ',outmode
       print*,'  Periodicity            : ',per
       print*,'  Time check             : ',trim(timecheck)
+      print*,'  Isentropic trajectories: ',trim(isen),thons
+      print*,'  Model-level trajs (2D) : ',trim(modlev)
       print*
+      if ( (isen.eq.'yes').and.(modlev.eq.'yes') ) then
+         print*,
+     >    '  WARNING: isentropic and 2D mode chosen -> 2D accepted'
+         print*
+         isen = 'no'
+      endif
+c     Init missing data value
+      mdv = -999.
       print*,'---- FIXED NUMERICAL PARAMETERS -------------------------'
       print*
       print*,'  Numerical scheme       : ',imethod
       print*,'  Number of iterations   : ',numit
       print*,'  Filename prefix        : ',prefix
       if (stat.ne.0) print*,'*** error allocating array wwt1 ***'
       allocate(p3t0(nx*ny*nz),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Pressure
       allocate(p3t1(nx*ny*nz),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
+      allocate(tht0(nx*ny*nz),stat=stat)
+      if (stat.ne.0) print*,'*** error allocating array tht0 ***'   ! Potential temperature
+      allocate(tht1(nx*ny*nz),stat=stat)
+      if (stat.ne.0) print*,'*** error allocating array tht1 ***'
+      allocate(sth0(nx*ny),stat=stat)
+      if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface potential temperature
+      allocate(sth1(nx*ny),stat=stat)
+      if (stat.ne.0) print*,'*** error allocating array sth1 ***'
 C     Get memory for trajectory arrays
       allocate(trainp(ntra,1,ncol),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array trainp   ***' ! Input start coordinates
       allocate(traout(ntra,ntim,4),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array yy0      ***' ! Y position (latitude)
       allocate(pp0(ntra),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array pp0      ***' ! Pressure
       allocate(leftflag(ntra),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array leftflag ***' ! Leaving-domain flag
+      allocate(theta(ntra),stat=stat)
+      if (stat.ne.0) print*,'*** error allocating array theta ***'    ! Potential temperature for isentropic trajectories
+      allocate(zindex(ntra),stat=stat)
+      if (stat.ne.0) print*,'*** error allocating array kind ***'     ! Vertical index for model-level trajectories
 c     Write some status information
       print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
       print*
       print*,'  xmin,xmax     : ',xmin,xmax
           enddo
          close(fid)
 c     Read start coordinates from trajectory file - check consistency of ref time
       else
-         print*,'Hallo', ntra,ncol
          call ropen_tra(cdfid,strname,ntra,1,ncol,reftmp,vars,inpmode)
          call read_tra (cdfid,trainp,ntra,1,ncol,inpmode)
          do i=1,ntra
             time   = trainp(i,1,1)
             xx0(i) = trainp(i,1,2)
       call hhmm2frac(tst,frac)
       tst = frac
       call hhmm2frac(ten,frac)
       ten = frac
-c     Check that all starting positions are above topography
+c     Get 3D and surface pressure from first data file
-      varname = 'P'
       filename = prefix//dat(1)
       call input_open (fid,filename)
+      if ( modlev.eq.'no' ) then
+         varname = 'P'
-      call input_grid
+         call input_grid
+     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
+     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
+      else
+         varname = 'P.ML'
+         call input_grid
-     >    (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
+     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
-     >     tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
+     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
+      endif
       call input_close(fid)
+c     Check that all starting positions are above topography
       do i=1,ntra
 C       Interpolate surface pressure to actual position (from first input file)
         x1 = xx0(i)
         y1 = yy0(i)
             leftflag(i) = 1
         endif
       enddo
+c     Special handling for isentropic trajectories - read potential
+c     temperature from S file or calculate it based on temperature and
+c     pressure from P file; then, calculate for each trajectory its
+c     potential temperature - which will stay fixed over time
+      if ( isen.eq.'yes' ) then
+c         Get potential temperature from S file
+          if ( thons.eq.1 ) then
+              filename = srefix//dat(1)
+              print*,' TH <- ',trim(filename)
+              call input_open (fid,filename)
+              varname='TH'
+              call input_wind
+     >            (fid,varname,tht1,tload,stagz,mdv,
+     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+              call input_close(fid)
+c         Calculate potential temperature from P file
+          else
+              filename = prefix//dat(1)
+              print*,' TH = T * (1000/P)^RDCP <- ',trim(filename)
+              call input_open (fid,filename)
+              varname='T'
+              call input_wind
+     >            (fid,varname,tht1,tload,stagz,mdv,
+     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+              call input_close(fid)
+              do i=1,nx*ny*nz
+                if (tht1(i).lt.100.) then
+                   tht1(i)=(tht1(i)+tzero)*( (1000./p3t1(i))**rdcp )
+                else
+                   tht1(i)=tht1(i)*( (1000./p3t1(i))**rdcp )
+                endif
+              enddo
+          endif
+c         Take surface potential temperature from lowest level
+          do i=1,nx*ny
+            sth1(i) = tht1(i)
+          enddo
+c         Calculate now the potential temperature of all trajectories
+          do i=1,ntra
+             x1 = xx0(i)
+             y1 = yy0(i)
+             p1 = pp0(i)
+             call get_index4 (xind,yind,pind,x1,y1,p1,0.,
+     >                 p3t1,p3t1,spt1,spt1,3,
+     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
+             theta(i) =
+     >        int_index4(tht1,tht1,nx,ny,nz,xind,yind,pind,0.,mdv)
+          enddo
+c         Write info about theta range of starting positions
+          thetamin = theta(1)
+          thetamax = theta(1)
+          do i=2,ntra
+            if ( theta(i).gt.thetamax ) thetamax = theta(i)
+            if ( theta(i).lt.thetamin ) thetamin = theta(i)
+          enddo
+c         Write some status information
+          print*
+          print*,
+     >     '---- THETA RANGE OF ISENTROPIC TRAJECTORIES -------------'
+          print*
+          print*,' Theta(min)             :',thetamin
+          print*,' Theta(max)             :',thetamax
+      endif
+c     Special handling for model-level (2D) trajectories - get the
+c     vertical index for each trajectory - which will remain fixed
+      if ( modlev.eq.'yes' ) then
+          do i=1,ntra
+             x1 = xx0(i)
+             y1 = yy0(i)
+             p1 = pp0(i)
+             call get_index4 (xind,yind,pind,x1,y1,p1,0.,
+     >                 p3t1,p3t1,spt1,spt1,3,
+     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
+             zindex(i) = pind
+          enddo
+          do i=1,nz
+            print*,i,p3t1(189+(141-1)*nx+(i-1)*nx*ny)
+          enddo
+             print*,x1,y1,p1
+             print*,xind,yind,pind
+c         Write info about zindex range of starting positions
+          zindexmin = zindex(1)
+          zindexmax = zindex(1)
+          do i=2,ntra
+            if ( zindex(i).gt.zindexmax ) zindexmax = zindex(i)
+            if ( zindex(i).lt.zindexmin ) zindexmin = zindex(i)
+          enddo
+c         Write some status information
+          print*
+          print*,
+     >     '---- INDEX RANGE OF MODEL-LEVEL TRAJECTORIES ------------'
+          print*
+          print*,' Zindex(min)            :',zindexmin
+          print*,' Zindex(max)            :',zindexmax
+      endif
 c     -----------------------------------------------------------------------
 c     Loop to calculate trajectories
 c     -----------------------------------------------------------------------
 c     Read wind fields and vertical grid from first file
       filename = prefix//dat(1)
       call frac2hhmm(tstart,tload)
-      write(*,'(a16,a20,f7.2)') '  (file,time) : ',
+      write(*,'(a16,a20,f9.2)') '  (file,time) : ',
      >                       trim(filename),tload
       call input_open (fid,filename)
       varname='U'                                      ! U
       call input_wind
      >    (fid,varname,uut1,tload,stagz,mdv,
      >     xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
       varname='V'                                      ! V
       call input_wind
      >    (fid,varname,vvt1,tload,stagz,mdv,
      >     xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+      if ( (modlev.eq.'no').and.(isen.eq.'no') ) then
-      varname='OMEGA'                                  ! OMEGA
+         varname='OMEGA'                                  ! OMEGA
-      call input_wind
+         call input_wind
-     >    (fid,varname,wwt1,tload,stagz,mdv,
+     >       (fid,varname,wwt1,tload,stagz,mdv,
-     >     xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+     >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+      endif
+      if ( modlev.eq.'no' ) then
+         call input_grid                                  ! GRID - AK,BK -> P
+     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
+     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
+      else
+         varname='P.ML'                                   ! GRID - P,PS
-      call input_grid                                  ! GRID
+         call input_grid                                  !
-     >    (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
+     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
-     >     tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
+     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
+      endif
       call input_close(fid)
+c     Special handling for isentropic trajectories - read potential
+c     temperature from S file or calculate it based on temperature and
+c     pressure from P file
+      if ( isen.eq.'yes' ) then
+c         Get potential temperature from S file
+          if ( thons.eq.1 ) then
+              filename = srefix//dat(1)
+              print*,' TH <- ',trim(filename)
+              call input_open (fid,filename)
+              varname='TH'
+              call input_wind
+     >            (fid,varname,tht1,tload,stagz,mdv,
+     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+              call input_close(fid)
+c         Calculate potential temperature from P file
+          else
+              filename = prefix//dat(1)
+              print*,' TH = T * (1000/P)^RDCP <- ',trim(filename)
+              call input_open (fid,filename)
+              varname='T'
+              call input_wind
+     >            (fid,varname,tht1,tload,stagz,mdv,
+     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+              call input_close(fid)
+              do i=1,nx*ny*nz
+                if (tht1(i).lt.100.) then
+                   tht1(i)=(tht1(i)+tzero)*( (1000./p3t1(i))**rdcp )
+                else
+                   tht1(i)=tht1(i)*( (1000./p3t1(i))**rdcp )
+                endif
+              enddo
+          endif
+c         Take surface potential temperature from lowest level
+          do i=1,nx*ny
+            sth1(i) = tht1(i)
+          enddo
+      endif
 c     Loop over all input files (time step is <timeinc>)
       do itm=1,numdat-1
 c       Calculate actual and next time
         time0 = tstart+real(itm-1)*timeinc*fbflag
         do i=1,nx*ny*nz
            uut0(i)=uut1(i)
            vvt0(i)=vvt1(i)
            wwt0(i)=wwt1(i)
            p3t0(i)=p3t1(i)
+           tht0(i)=tht1(i)
         enddo
         do i=1,nx*ny
            spt0(i)=spt1(i)
+           sth0(i)=sth1(i)
         enddo
 c       Read wind fields and surface pressure at next time
         filename = prefix//dat(itm+1)
         call frac2hhmm(time1,tload)
-        write(*,'(a16,a20,f7.2)') '  (file,time) : ',
+        write(*,'(a16,a20,f9.2)') '  (file,time) : ',
      >                          trim(filename),tload
         call input_open (fid,filename)
         varname='U'                                     ! U
         call input_wind
      >       (fid,varname,uut1,tload,stagz,mdv,
      >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
         varname='V'                                     ! V
         call input_wind
      >       (fid,varname,vvt1,tload,stagz,mdv,
      >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+        if ( (modlev.eq.'no').and.(isen.eq.'no') ) then
-        varname='OMEGA'                                ! OMEGA
+           varname='OMEGA'                              ! OMEGA
-        call input_wind
+           call input_wind
-     >       (fid,varname,wwt1,tload,stagz,mdv,
+     >          (fid,varname,wwt1,tload,stagz,mdv,
-     >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+     >           xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+        endif
+        if ( modlev.eq.'no' ) then
+           call input_grid                                  ! GRID - AK,NK -> P
+     >          (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
+     >           tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
+        else
+           varname='P.ML'                                   ! GRID - P,PS
-        call input_grid                                ! GRID
+           call input_grid                                  !
      >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
      >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
+        endif
         call input_close(fid)
+c     Special handling for isentropic trajectories - read potential
+c     temperature from S file or calculate it based on temperature and
+c     pressure from P file
+      if ( isen.eq.'yes' ) then
+c         Get TH from S file
+          if ( thons.eq.1 ) then
+              filename = srefix//dat(itm+1)
+              print*,' TH <- ',trim(filename)
+              call input_open (fid,filename)
+              varname='TH'
+              call input_wind
+     >            (fid,varname,tht1,tload,stagz,mdv,
+     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+              call input_close(fid)
+c         Calculate potential temperature from P file
+          else
+              filename = prefix//dat(itm+1)
+              print*,' TH = T * (1000/P)^RDCP <- ',trim(filename)
+              call input_open (fid,filename)
+              varname='T'
+              call input_wind
+     >            (fid,varname,tht1,tload,stagz,mdv,
+     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
+              call input_close(fid)
+              do i=1,nx*ny*nz
+                if (tht1(i).lt.100.) then
+                   tht1(i)=(tht1(i)+tzero)*( (1000./p3t1(i))**rdcp )
+                else
+                   tht1(i)=tht1(i)*( (1000./p3t1(i))**rdcp )
+                endif
+              enddo
+          endif
+c         Take surface potential temperature from lowest level
+          do i=1,nx*ny
+            sth1(i) = tht1(i)
+          enddo
+      endif
 C       Determine the first and last loop indices
         if (numdat.eq.2) then
           filo = nint(tst/ts)+1
           lalo = nint((timeinc-ten)/ts)
           do i=1,ntra
 C           Check if trajectory has already left the data domain
             if (leftflag(i).ne.1) then
-c             Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
+c             3D: Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
-              if (imethod.eq.1) then
+              if ( (imethod.eq.1   ).and.
+     >             (isen   .eq.'no').and.
+     >             (modlev .eq.'no') )
+     >        then
-                 call euler(
+                 call euler_3d(
      >               xx1,yy1,pp1,leftflag(i),
      >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
      >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
      >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
      >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
-c             Runge-Kutta timestep (x0,y0,p0 -> x1,y1,p1)
+c             3D: Runge-Kutta timestep (x0,y0,p0 -> x1,y1,p1)
-              else if (imethod.eq.2) then
+              else if ( (imethod.eq.2   ).and.
+     >                  (isen   .eq.'no').and.
+     >                  (modlev .eq.'no') )
+     >        then
                  call runge(
      >               xx1,yy1,pp1,leftflag(i),
      >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
      >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
      >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
      >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
+c             ISENTROPIC: Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
+              else if ( (imethod.eq.1    ).and.
+     >                  (isen   .eq.'yes').and.
+     >                  (modlev .eq.'no' ) )
+     >        then
+                 call euler_isen(
+     >               xx1,yy1,pp1,leftflag(i),
+     >               xx0(i),yy0(i),pp0(i),theta(i),reltpos0,reltpos1,
+     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
+     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,
+     >               sth0,sth1,tht0,tht1,
+     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
+c             MODEL-LEVEL (2D): Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
+              else if ( (imethod.eq.1    ).and.
+     >                  (isen   .eq.'no' ).and.
+     >                  (modlev .eq.'yes') )
+     >        then
+                 call euler_2d(
+     >               xx1,yy1,pp1,leftflag(i),
+     >               xx0(i),yy0(i),pp0(i),zindex(i),reltpos0,reltpos1,
+     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
+     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,
+     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
               endif
 c             Update trajectory position, or increase number of trajectories leaving domain
               if (leftflag(i).eq.1) then
                 leftcount=leftcount+1
         enddo
       enddo
-      print*,(traout(1,1,i),i=1,4)
 c     Write trajectory file
       vars(1)  ='time'
       vars(2)  ='lon'
       vars(3)  ='lat'
       vars(4)  ='p'
 C     numit	 input	number of iterations
 C     jump	 input  flag (=1 trajectories don't enter the ground)
 C     left	 output	flag (=1 if trajectory leaves data domain)
 c     -------------------------------------------------------------------
-c     Iterative Euler time step
+c     Iterative Euler time step (KINEMATIC 3D TRAJECTORIES)
 c     -------------------------------------------------------------------
-      subroutine euler(x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
+      subroutine euler_3d(x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
-     >                 deltat,numit,jump,mdv,wfactor,fbflag,
+     >                deltat,numit,jump,mdv,wfactor,fbflag,
-     >		       spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
+     >		          spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
-     >                 xmin,ymin,dx,dy,per,hem,nx,ny,nz)
+     >                xmin,ymin,dx,dy,per,hem,nx,ny,nz)
       implicit none
 c     Declaration of subroutine parameters
       integer      nx,ny,nz
 C       Check if trajectory leaves data domain
         if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
      >       ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
      >         (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
      >  then
+          left=1
+          goto 100
+        endif
+      enddo
+c     Exit point for subroutine
+continue
+      return
+      end
+c     -------------------------------------------------------------------
+c     Iterative Euler time step (ISENTROPIC)
+c     -------------------------------------------------------------------
+      subroutine euler_isen(x1,y1,p1,left,x0,y0,p0,theta,
+     >                reltpos0,reltpos1,
+     >                deltat,numit,jump,mdv,wfactor,fbflag,
+     >                spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
+     >                sth0,sth1,tht0,tht1,
+     >                xmin,ymin,dx,dy,per,hem,nx,ny,nz)
+      implicit none
+c     Declaration of subroutine parameters
+      integer      nx,ny,nz
+      real         x1,y1,p1
+      integer      left
+      real         x0,y0,p0
+      real         reltpos0,reltpos1
+      real         deltat
+      integer      numit
+      integer      jump
+      real         wfactor
+      integer      fbflag
+      real         spt0(nx*ny)   ,spt1(nx*ny)
+      real         sth0(nx*ny)   ,sth1(nx*ny)
+      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
+      real         vvt0(nx*ny*nz),vvt1(nx*ny*nz)
+      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
+      real         tht0(nx*ny*nz),tht1(nx*ny*nz)
+      real         xmin,ymin,dx,dy
+      real         per
+      integer      hem
+      real         mdv
+      real         theta
+c     Numerical and physical constants
+      real         deltay
+      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
+      real         pi
+      parameter    (pi=3.1415927)    ! Pi
+c     Auxiliary variables
+      real         xmax,ymax
+      real         xind,yind,pind
+      real         u0,v0,w0,u1,v1,w1,u,v,w,sp
+      integer      icount
+      character    ch
+c     Externals
+      real         int_index4
+      external     int_index4
+c     Reset the flag for domain-leaving
+      left=0
+c     Set the esat-north bounray of the domain
+      xmax = xmin+real(nx-1)*dx
+      ymax = ymin+real(ny-1)*dy
+C     Interpolate wind fields to starting position (x0,y0,p0)
+      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
+     >                 p3d0,p3d1,spt0,spt1,3,
+     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
+      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
+      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
+C     For first iteration take ending position equal to starting position
+      x1=x0
+      y1=y0
+      p1=p0
+C     Iterative calculation of new position
+      do icount=1,numit
+C        Calculate new winds for advection
+         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
+     >                    p3d0,p3d1,spt0,spt1,3,
+     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
+         u1 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
+         v1 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
+c        Get the new velocity in between
+         u=(u0+u1)/2.
+         v=(v0+v1)/2.
+C        Calculate new positions
+         x1 = x0 + fbflag*u*deltat/(deltay*cos(y0*pi/180.))
+         y1 = y0 + fbflag*v*deltat/deltay
+c        Get the pressure on the isentropic surface at the new position
+         call get_index4 (xind,yind,pind,x1,y1,theta,reltpos1,
+     >                    tht0,tht1,sth0,sth1,1,
+     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
+         p1 = int_index4(p3d0,p3d1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
+c       Handle pole problems (crossing and near pole trajectory)
+        if ((hem.eq.1).and.(y1.gt.90.)) then
+          y1=180.-y1
+          x1=x1+per/2.
+        endif
+        if ((hem.eq.1).and.(y1.lt.-90.)) then
+          y1=-180.-y1
+          x1=x1+per/2.
+        endif
+        if (y1.gt.89.99) then
+           y1=89.99
+        endif
+c       Handle crossings of the dateline
+        if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
+           x1=xmin+amod(x1-xmin,per)
+        endif
+        if ((hem.eq.1).and.(x1.lt.xmin)) then
+           x1=xmin+per+amod(x1-xmin,per)
+        endif
+C       Interpolate surface pressure to actual position
+        call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
+     >                   p3d0,p3d1,spt0,spt1,3,
+     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
+        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+c       Handle trajectories which cross the lower boundary (jump flag)
+        if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
+C       Check if trajectory leaves data domain
+        if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
+     >       ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
+     >         (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
+     >  then
+          left=1
+          goto 100
+        endif
+      enddo
+c     Exit point for subroutine
+continue
+      return
+      end
+c     -------------------------------------------------------------------
+c     Iterative Euler time step (MODEL-LEVEL, 2D)
+c     -------------------------------------------------------------------
+      subroutine euler_2d(x1,y1,p1,left,x0,y0,p0,zindex,
+     >                reltpos0,reltpos1,
+     >                deltat,numit,jump,mdv,wfactor,fbflag,
+     >                spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
+     >                xmin,ymin,dx,dy,per,hem,nx,ny,nz)
+      implicit none
+c     Declaration of subroutine parameters
+      integer      nx,ny,nz
+      real         x1,y1,p1
+      integer      left
+      real         x0,y0,p0
+      real         reltpos0,reltpos1
+      real         deltat
+      integer      numit
+      integer      jump
+      real         wfactor
+      integer      fbflag
+      real         spt0(nx*ny)   ,spt1(nx*ny)
+      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
+      real         vvt0(nx*ny*nz),vvt1(nx*ny*nz)
+      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
+      real         xmin,ymin,dx,dy
+      real         per
+      integer      hem
+      real         mdv
+      real         zindex
+c     Numerical and physical constants
+      real         deltay
+      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
+      real         pi
+      parameter    (pi=3.1415927)    ! Pi
+      real         eps
+      parameter    (eps=0.001)
+c     Auxiliary variables
+      real         xmax,ymax
+      real         xind,yind,pind
+      real         u0,v0,w0,u1,v1,w1,u,v,w,sp
+      integer      icount
+      character    ch
+c     Externals
+      real         int_index4
+      external     int_index4
+c     Reset the flag for domain-leaving
+      left=0
+c     Set the esat-north bounray of the domain
+      xmax = xmin+real(nx-1)*dx
+      ymax = ymin+real(ny-1)*dy
+C     Interpolate wind fields to starting position (x0,y0,p0)
+      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
+     >                 p3d0,p3d1,spt0,spt1,3,
+     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
+      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,zindex,reltpos0,mdv)
+      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,zindex,reltpos0,mdv)
+C     For first iteration take ending position equal to starting position
+      x1=x0
+      y1=y0
+      p1=p0
+C     Iterative calculation of new position
+      do icount=1,numit
+C        Calculate new winds for advection
+         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
+     >                    p3d0,p3d1,spt0,spt1,3,
+     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
+         u1=int_index4(uut0,uut1,nx,ny,nz,xind,yind,zindex,reltpos1,mdv)
+         v1=int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,zindex,reltpos1,mdv)
+         if ( abs(u1-mdv).lt.eps ) then
+             left = 1
+             goto 100
+         endif
+         if ( abs(v1-mdv).lt.eps ) then
+             left = 1
+             goto 100
+         endif
+c        Get the new velocity in between
+         u=(u0+u1)/2.
+         v=(v0+v1)/2.
+C        Calculate new positions
+         x1 = x0 + fbflag*u*deltat/(deltay*cos(y0*pi/180.))
+         y1 = y0 + fbflag*v*deltat/deltay
+c        Get the pressure on the model surface at the new position
+         xind = (x1 - xmin ) / dx + 1.
+         yind = (y1 - ymin ) / dy + 1.
+         p1 =
+     >     int_index4(p3d0,p3d1,nx,ny,nz,xind,yind,zindex,reltpos1,mdv)
+         if ( abs(p1-mdv).lt.eps ) then
+             left = 1
+             goto 100
+         endif
+c       Handle pole problems (crossing and near pole trajectory)
+        if ((hem.eq.1).and.(y1.gt.90.)) then
+          y1=180.-y1
+          x1=x1+per/2.
+        endif
+        if ((hem.eq.1).and.(y1.lt.-90.)) then
+          y1=-180.-y1
+          x1=x1+per/2.
+        endif
+        if (y1.gt.89.99) then
+           y1=89.99
+        endif
+c       Handle crossings of the dateline
+        if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
+           x1=xmin+amod(x1-xmin,per)
+        endif
+        if ((hem.eq.1).and.(x1.lt.xmin)) then
+           x1=xmin+per+amod(x1-xmin,per)
+        endif
+C       Interpolate surface pressure to actual position
+        call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
+     >                   p3d0,p3d1,spt0,spt1,3,
+     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
+        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+c       Handle trajectories which cross the lower boundary (jump flag)
+        if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
+C       Check if trajectory leaves data domain
+        if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
+     >       ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
+     >         (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
+     >  then
           left=1
           goto 100
         endif
       enddo

Subversion Repositories lagranto.ecmwf

(root)/trunk/caltra/caltra.f @ 33 – Rev 13 → 19