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

       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*10		              	     gridflag        ! "ideal" or "real"
+      real                                   aglconst        ! constant level above ground
+      real                                   isen            ! isentropic level
+      character*80                           vertmode        ! Vertical mode (kinematic,aglconst)
 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 (:)     :: 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       ! 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                                wstep
       character*80                           radunit
       real                                   lon1,lat1,lon2,lat2
       real                                   scale_max,scale_min
+c     Externals
+      real         int_index4
+      external     int_index4
 c     --------------------------------------------------------------------
 c     Start of program, Read parameters
 c     --------------------------------------------------------------------
 c     Write start message
       read(9,*) reftime(6)                  ! time range (in min)
 c     Read flag for 'no time check'
       read(9,*) timecheck
+c     Read flag for trajectories at constant level above ground
+      read(9,*) aglconst
+c     Read flag for trajectories at constant isentropic level
+      read(9,*) isen
 c     Close the input file
       close(9)
+c     Set the vertical mode of the trajectory calculation
+      vertmode = 'kinematic'
+      if ( aglconst.gt.0. ) then
+         vertmode = 'aglconst'
+      endif
+      if ( isen.gt.0. ) then
+         vertmode = 'isen'
+      endif
 c     Calculate the number of output time steps
       ntim = abs(reftime(6)/deltout) + 1
 c     Set the formats of the input and output files
       call mode_tra(inpmode,strname)
          print*,'  Input format           : ',inpmode
       endif
       print*,'  Output format          : ',outmode
       print*,'  Periodicity            : ',per
       print*,'  Time check             : ',trim(timecheck)
+      if ( vertmode.eq.'aglconst' ) then
+         print*,'  Z = const AGL          : ',aglconst
+      endif
+      if ( vertmode.eq.'isen' ) then
+         print*,'  isentropic             : ',isen
+      endif
       print*
       print*,'---- FIXED NUMERICAL PARAMETERS -------------------------'
       print*
       print*,'  Number of iterations   : ',numit
       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
+c     Get memory for potential temperature (isentropic trajectories)
+      if ( vertmode.eq.'isen' ) then
+          allocate(tht0(nx*ny*nz),stat=stat)
+          if (stat.ne.0) print*,'*** error allocating array tht0 ***'
+          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 sth0 ***'
+          allocate(sth1(nx*ny),stat=stat)
+          if (stat.ne.0) print*,'*** error allocating array sth1 ***'
+      endif
 c     Write some status information
       print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
       print*
       print*,'  xmin,xmax     : ',xmin,xmax
       print*,'  ymin,ymax     : ',ymin,ymax
       print*,'                (hour)  :',reftime(4)
       print*,'                (min)   :',reftime(5)
       print*,' Time range             :',reftime(6),' min'
       print*
-C     Save starting positions
+C     Save starting positions - vertical position will be corrected below for <aglconst> and <isen>
       itim = 1
       do i=1,ntra
          traout(i,itim,1) = 0.
          traout(i,itim,2) = xx0(i)
          traout(i,itim,3) = yy0(i)
       call hhmm2frac(tst,frac)
       tst = frac
       call hhmm2frac(ten,frac)
       ten = frac
+c     Adjust the starting heights in mode <aglconst>
+      if ( vertmode.eq.'aglconst' ) then
+        filename = prefix//dat(1)
+        call input_open (fid,filename)
+	      varname='geopot'				! geopot.height
+	      call input_var_wrf(fid,varname,p3t0,nx,ny,nz)
+	      varname='geopots'				! geopot.height at surface
+	      call input_var_wrf(fid,varname,spt0,nx,ny,1)
+        call input_close(fid)
+        filename = prefix//dat(2)
+        call input_open (fid,filename)
+	      varname='geopot'				! geopot.height
+	      call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
+	      varname='geopots'				! geopot.height at surface
+	      call input_var_wrf(fid,varname,spt1,nx,ny,1)
+        call input_close(fid)
+        do i=1,ntra
+           call get_index4 (xind,yind,pind,xx0(i),yy0(i),0.,reltpos0,
+     >                      p3t0,p3t1,spt0,spt1,3,
+     >                      nx,ny,nz,xmin,ymin,dx,dy,mdv)
+           sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos0,mdv)
+           traout(i,1,4) = sp + aglconst
+        enddo
+      endif
+c     Adjust the starting heights in mode <isen>
+      if ( vertmode.eq.'isen' ) then
+        filename = prefix//dat(1)
+        call input_open (fid,filename)
+	      varname='geopot'				! geopot.height
+	      call input_var_wrf(fid,varname,p3t0,nx,ny,nz)
+	      varname='geopots'				! geopot.height at surface
+	      call input_var_wrf(fid,varname,spt0,nx,ny,1)
+          varname='TH'					! TH
+          call input_var_wrf(fid,varname,tht0,nx,ny,nz)
+          do i=1,nx*ny
+              sth0(i) = tht0(i)
+          enddo
+        call input_close(fid)
+        filename = prefix//dat(2)
+        call input_open (fid,filename)
+	      varname='geopot'				! geopot.height
+	      call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
+	      varname='geopots'				! geopot.height at surface
+	      call input_var_wrf(fid,varname,spt1,nx,ny,1)
+          varname='TH'					! TH
+          call input_var_wrf(fid,varname,tht1,nx,ny,nz)
+          do i=1,nx*ny
+              sth1(i) = tht1(i)
+          enddo
+        call input_close(fid)
+        do i=1,ntra
+          call get_index4 (xind,yind,pind,xx0(i),yy0(i),isen,reltpos0,
+     >                    tht0,tht1,sth0,sth1,1,
+     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
+          traout(i,1,4) = int_index4
+     >          (p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
+        enddo
+      endif
 c     -----------------------------------------------------------------------
 c     Loop to calculate trajectories
 c     -----------------------------------------------------------------------
 c     Write some status information
 	  call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
 	  varname='geopots'				! geopot.height at surface
 	  call input_var_wrf(fid,varname,spt1,nx,ny,1)
+      if (vertmode.eq.'isen') then
+           varname='TH'					! TH
+ 	       call input_var_wrf(fid,varname,tht1,nx,ny,nz)
+           do i=1,nx*ny
+              sth1(i) = tht1(i)
+           enddo
+      endif
       call input_close(fid)
 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
         time1 = time0+timeinc*fbflag
 c       Copy old velocities and pressure fields to new ones
         do i=1,nx*ny*nz
            uut0(i)=uut1(i)
            vvt0(i)=vvt1(i)
            wwt0(i)=wwt1(i)
         enddo
         do i=1,nx*ny
            spt0(i)=spt1(i)
         enddo
+c       Copy potential temperature
+        if ( vertmode.eq.'isen' ) then
+           do i=1,nx*ny*nz
+              tht0(i)=tht1(i)
+           enddo
+           do i=1,nx*ny
+             sth0(i)=sth1(i)
+           enddo
+        endif
 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) : ',
 	    call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
 	    varname='geopots'				! geopot.height
 	    call input_var_wrf(fid,varname,spt1,nx,ny,1)
+        if (vertmode.eq.'isen') then
+           varname='TH'					! TH
+ 	       call input_var_wrf(fid,varname,tht1,nx,ny,nz)
+           do i=1,nx*ny
+              sth1(i) = tht1(i)
+           enddo
+        endif
         call input_close(fid)
 C       Determine the first and last loop indices
         if (numdat.eq.2) then
           filo = nint(tst/ts)+1
 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)
+              if ( vertmode.eq.'kinematic' ) then
-              call euler(
+                 call euler_kinematic(
      >               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,mpx,mpy)
+              elseif ( vertmode.eq.'aglconst' ) then
+                 call euler_aglconst(
+     >               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,
+     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,aglconst)
+              elseif ( vertmode.eq.'isen' ) then
+                 call euler_isen(
+     >               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,
+     >               tht0,tht1,sth0,sth1,
+     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,isen)
+              endif
 c             Update trajectory position, or increase number of trajectories leaving domain
               if (leftflag(i).eq.1) then
                 leftcount=leftcount+1
                 if ( leftcount.lt.10 ) then
 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 trajectory
 c     -------------------------------------------------------------------
+      subroutine euler_kinematic
-      subroutine euler(x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
+     >          (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,mpx,mpy)
+     >           xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy)
       implicit none
       return
       end
+c     -------------------------------------------------------------------
+c     Iterative Euler time step - constant level above ground
+c     -------------------------------------------------------------------
+      subroutine euler_aglconst
+     >          (x1,y1,p1,left,x0,y0,p0,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,mpx,mpy,aglconst)
+      implicit none
+c     Flag for test mode
+      integer      test
+      parameter    (test=0)
+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
+      integer      per
+      integer      hem
+      real         mdv
+      real         mpx(nx*ny),mpy(nx*ny)
+      real         aglconst
+c     Auxiliary variables
+      real         xmax,ymax
+      real	       xind,yind,pind
+      real	       u0,v0,u1,v1,u,v,w,sp
+      integer	   icount
+      character    ch
+      real         mpsc_x,mpsc_y
+c     Externals
+      real         int_index4
+      external     int_index4
+c     Reset the flag for domain-leaving
+      left=0
+c     Set the esat-north boundary of the domain
+      xmax = xmin+real(nx-1)*dx
+      ymax = ymin+real(ny-1)*dy
+C     Set parcel height to aglconst above topography
+      call get_index4 (xind,yind,pind,x1,y1,0.,reltpos0,
+     >                 p3d0,p3d1,spt0,spt1,3,
+     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
+      sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos0,mdv)
+      p0 = sp + aglconst
+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        Get the mapping scale factors for this position
+         mpsc_x = int_index4 (mpx,mpx,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+         mpsc_y = int_index4 (mpy,mpy,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+C        Calculate new positions (adapted for cartesian grid)
+         x1 = x0 + fbflag * deltat * u * dx/mpsc_x
+         y1 = y0 + fbflag * deltat * v * dy/mpsc_y
+c        Write the update positions for tests
+         if ( (icount.eq.3).and.(test.eq.1) ) then
+            write(*,'(10f10.2)') x0,u,x1-x0,p0,w,p1-p0
+         endif
+C       Set trajectory height to aglconst above topography
+        call get_index4 (xind,yind,pind,x1,y1,0.,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)
+        p1 = sp + aglconst
+c       Handle trajectories which cross the lower boundary (jump flag)
+        if ((jump.eq.1).and.(p1.lt.sp)) then
+            print*,'jump'
+            p1=sp+10.
+        endif
+c       Handle periodioc boundaries for 'ideal' mode
+        if ( per.eq.1 ) then
+	        if ( x1.gt.xmax ) x1=x1-xmax
+	        if ( x1.lt.xmin ) x1=x1+xmax
+        endif
+C       Check if trajectory leaves data domain
+        if ( (x1.lt.xmin).or.(x1.gt.xmax).or.
+     >       (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.lt.sp) )   ! geopot : .lt. ; pressure: .gt.
+     >  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,reltpos0,reltpos1,
+     >           deltat,numit,jump,mdv,wfactor,fbflag,
+     >		     spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
+     >           tht0,tht1,sth0,sth1,
+     >           xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,isen)
+      implicit none
+c     Flag for test mode
+      integer      test
+      parameter    (test=0)
+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         tht0(nx*ny*nz),tht1(nx*ny*nz),sth0(nx*ny),sth1(nx*ny)
+      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
+      real         xmin,ymin,dx,dy
+      integer      per
+      integer      hem
+      real         mdv
+      real         mpx(nx*ny),mpy(nx*ny)
+      real         isen
+c     Auxiliary variables
+      real         xmax,ymax
+      real	       xind,yind,pind
+      real	       u0,v0,u1,v1,u,v,w,sp
+      integer	   icount
+      character    ch
+      real         mpsc_x,mpsc_y
+c     Externals
+      real         int_index4
+      external     int_index4
+c     Reset the flag for domain-leaving
+      left=0
+c     Set the esat-north boundary 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        Get the mapping scale factors for this position
+         mpsc_x = int_index4 (mpx,mpx,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+         mpsc_y = int_index4 (mpy,mpy,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+C        Calculate new positions (adapted for cartesian grid)
+         x1 = x0 + fbflag * deltat * u * dx/mpsc_x
+         y1 = y0 + fbflag * deltat * v * dy/mpsc_y
+c        Write the update positions for tests
+         if ( (icount.eq.3).and.(test.eq.1) ) then
+            write(*,'(10f10.2)') x0,u,x1-x0,p0,w,p1-p0
+         endif
+c        Set trajectory height to isentropic level
+         call get_index4 (xind,yind,pind,x1,y1,isen,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       Interpolate surface geop.height to actual position
+        call get_index4 (xind,yind,pind,x1,y1,0.,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.lt.sp)) then
+            print*,'jump'
+            p1=sp+10.
+        endif
+c       Handle periodioc boundaries for 'ideal' mode
+        if ( per.eq.1 ) then
+	        if ( x1.gt.xmax ) x1=x1-xmax
+	        if ( x1.lt.xmin ) x1=x1+xmax
+        endif
+C       Check if trajectory leaves data domain
+        if ( (x1.lt.xmin).or.(x1.gt.xmax).or.
+     >       (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.lt.sp) )   ! geopot : .lt. ; pressure: .gt.
+     >  then
+          left=1
+          goto 100
+        endif
+      enddo
+c     Exit point for subroutine
+continue
+      return
+      end
 c     ----------------------------------------------------------------------
 c     Get spherical distance between lat/lon points
 c     ----------------------------------------------------------------------

Subversion Repositories lagranto.wrf

(root)/trunk/caltra/caltra.f – Rev 19 → 21