Subversion Repositories lagranto.ecmwf

Rev

Rev 15 | Rev 44 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 15 Rev 39
Line 63... Line 63...
63 
  character(len=80) :: tvar(200)       ! Tracing variable name (only the variable)
 63 
  character(len=80) :: tvar(200)       ! Tracing variable name (only the variable)
64 
  character(len=80) :: tfil(200)       ! Filename prefix
 64 
  character(len=80) :: tfil(200)       ! Filename prefix
65 
  real :: fac(200)                     ! Scaling factor
 65 
  real :: fac(200)                     ! Scaling factor
66 
  real :: shift_val(200)               ! Shift in space and time relative to trajectory position
 66 
  real :: shift_val(200)               ! Shift in space and time relative to trajectory position
67 
  character(len=80) :: shift_dir(200)  ! Direction of shift
 67 
  character(len=80) :: shift_dir(200)  ! Direction of shift
- 
 
 68 
  character(len=80) :: shift_rel(200)  ! Operator/relator for variable
68 
  integer :: compfl(200)               ! Computation flag (1=compute)
 69 
  integer :: compfl(200)               ! Computation flag (1=compute)
69 
  integer :: numdat                    ! Number of input files
 70 
  integer :: numdat                    ! Number of input files
70 
  character(len=11) :: dat(ndatmax)    ! Dates of input files
 71 
  character(len=11) :: dat(ndatmax)    ! Dates of input files
71 
  real :: timeinc                      ! Time increment between input files
 72 
  real :: timeinc                      ! Time increment between input files
72 
  real :: tst                          ! Time shift of start relative to first data file
 73 
  real :: tst                          ! Time shift of start relative to first data file
Line 759... Line 760...
759 
  endif
 760 
  endif
760 
  ! Bojan
 761 
  ! Bojan
761 
 
 762 
 
762 
  ! Split the tracing variables
 763 
  ! Split the tracing variables
763 
  do i=1,ntrace0
 764 
  do i=1,ntrace0
764 
    call splitvar(tvar(i),shift_val(i),shift_dir(i) )
 765 
    call splitvar(tvar(i),shift_val(i),shift_dir(i),shift_rel(i) )
765 
  enddo
 766 
  enddo
766 
 
 767 
 
767 
 
 768 
 
768 
  ! Split the variable name and set flags
 769 
  ! Split the variable name and set flags
769 
  do i=ntrace0+1,ntrace1
 770 
  do i=ntrace0+1,ntrace1
770 
 
 771 
 
771 
     ! Set the scaling factor
 772 
     ! Set the scaling factor
772 
     fac(i) = 1.
 773 
     fac(i) = 1.
773 
 
 774 
 
774 
     ! Set the base variable name, the shift and the direction
 775 
     ! Set the base variable name, the shift and the direction
775 
     call splitvar(tvar(i),shift_val(i),shift_dir(i) )
 776 
     call splitvar(tvar(i),shift_val(i),shift_dir(i),shift_rel(i) )
776 
 
 777 
 
777 
     ! Set the prefix of the file name for additional fields
 778 
     ! Set the prefix of the file name for additional fields
778 
     tfil(i)='*'
 779 
     tfil(i)='*'
779 
     do j=1,n_pvars
 780 
     do j=1,n_pvars
780 
        if ( tvar(i).eq.pvars(j) ) tfil(i)=charp
 781 
        if ( tvar(i).eq.pvars(j) ) tfil(i)=charp
Line 807... Line 808...
807 
          ( shift_dir(i).ne.'HPA(ABS)').and. &
 808 
          ( shift_dir(i).ne.'HPA(ABS)').and. &
808 
          ( shift_dir(i).ne.'KM(LON)' ).and. &
 809 
          ( shift_dir(i).ne.'KM(LON)' ).and. &
809 
          ( shift_dir(i).ne.'KM(LAT)' ).and. &
 810 
          ( shift_dir(i).ne.'KM(LAT)' ).and. &
810 
          ( shift_dir(i).ne.'H'       ).and. &
 811 
          ( shift_dir(i).ne.'H'       ).and. &
811 
          ( shift_dir(i).ne.'MIN'     ).and. &
 812 
          ( shift_dir(i).ne.'MIN'     ).and. &
- 
 
 813 
          ( shift_dir(i).ne.'INDP'    ).and. &
- 
 
 814 
          ( shift_dir(i).ne.'PMIN'    ).and. &
812 
          ( shift_dir(i).ne.'INDP'    ) ) then
 815 
          ( shift_dir(i).ne.'PMAX'    ) ) then
813 
        print*,' ERROR: shift mode ',trim(shift_dir(i)), ' not supported'
 816 
        print*,' ERROR: shift mode ',trim(shift_dir(i)), ' not supported'
814 
        stop
 817 
        stop
815 
     endif
 818 
     endif
816 
  enddo
 819 
  enddo
817 
 
 820 
 
Line 819... Line 822...
819 
  print*
 822 
  print*
820 
  print*,'---- COMPLETE TABLE FOR TRACING -------------------------'
 823 
  print*,'---- COMPLETE TABLE FOR TRACING -------------------------'
821 
  print*
 824 
  print*
822 
  do i=1,ntrace1
 825 
  do i=1,ntrace1
823 
     if ( ( shift_dir(i).ne.'nil' ) ) then
 826 
     if ( ( shift_dir(i).ne.'nil' ) ) then
824 
        write(*,'(i4,a4,a8,f10.2,a8,3x,a4,i5)') i,' : ',trim(tvar(i)), &
 827 
        write(*,'(i4,a4,a8,f10.2,a8,a8,3x,a4,i5)') i,' : ',trim(tvar(i)), &
825 
             shift_val(i),trim(shift_dir(i)),tfil(i),compfl(i)
 828 
             shift_val(i),trim(shift_dir(i)),trim(shift_rel(i)),tfil(i),compfl(i)
826 
     else
 829 
     else
827 
        write(*,'(i4,a4,a8,10x,8x,3x,a4,i5)') &
 830 
        write(*,'(i4,a4,a8,10x,8x,3x,a4,i5)') &
828 
             i,' : ',trim(tvar(i)),tfil(i),compfl(i)
 831 
             i,' : ',trim(tvar(i)),tfil(i),compfl(i)
829 
     endif
 832 
     endif
830 
  enddo
 833 
  enddo
Line 929... Line 932...
929 
            filename = trim(tfil(i))//trim(dat(itime0))
 932 
            filename = trim(tfil(i))//trim(dat(itime0))
930 
            call frac2hhmm(time0,tload)
 933 
            call frac2hhmm(time0,tload)
931 
            varname  = tvar(i)
 934 
            varname  = tvar(i)
932 
            write(*,'(a23,a20,a3,a5,f7.2)') '    ->  loading          : ', trim(filename),'  ',trim(varname),tload
 935 
            write(*,'(a23,a20,a3,a5,f7.2)') '    ->  loading          : ', trim(filename),'  ',trim(varname),tload
933 
            call input_open (fid,filename)
 936 
            call input_open (fid,filename)
- 
 
 937 
            mdv = -999.
934 
            call input_wind &
 938 
            call input_wind &
935 
                 (fid,varname,f3t0,tload,stagz,mdv, &
 939 
                 (fid,varname,f3t0,tload,stagz,mdv, &
936 
                 xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
 940 
                 xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
937 
            call input_grid &
 941 
            call input_grid &
938 
                 (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny, &
 942 
                 (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny, &
Line 1057... Line 1061...
1057 
            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
 1061 
            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
1058 
 
 1062 
 
1059 
               if ( pind.gt.nz ) then ! pressure too low, index too high
 1063 
               if ( pind.gt.nz ) then ! pressure too low, index too high
1060 
                 err_c1 = err_c1 + 1
 1064 
                 err_c1 = err_c1 + 1
1061 
                 if ( err_c1.lt.10 ) then
 1065 
                 if ( err_c1.lt.10 ) then
1062 
                    write(*,'(Af5.3A)') ' WARNING: pressure too low (pind = ',pind,') => adjusted to highest level (pind=nz.)'
 1066 
                    write(*,'(a,f7.3,a)') ' WARNING: pressure too low (pind = ',pind,') => adjusted to highest level (pind=nz.)'
1063 
                    print*,'(x0,y0,p0)=',x0,y0,p0
 1067 
                    print*,'(x0,y0,p0)=',x0,y0,p0
1064 
                    pind = real(nz)
 1068 
                    pind = real(nz)
1065 
                 elseif ( err_c1.eq.10 ) then
 1069 
                 elseif ( err_c1.eq.10 ) then
1066 
                    print*,' WARNING: more pressures too low -> adjusted to highest level '
 1070 
                    print*,' WARNING: more pressures too low -> adjusted to highest level '
1067 
                    pind = real(nz)
 1071 
                    pind = real(nz)
Line 1070... Line 1074...
1070 
                 endif
 1074 
                 endif
1071
1075
1072 
               elseif (pind.lt.1.) then ! pressure too high, index too low
 1076 
               elseif (pind.lt.1.) then ! pressure too high, index too low
1073 
                 err_c2 = err_c2 + 1
 1077 
                 err_c2 = err_c2 + 1
1074 
                 if ( err_c2.lt.10 ) then
 1078 
                 if ( err_c2.lt.10 ) then
1075 
                    write(*,'(Af5.3A)') ' WARNING: pressure too high (pind = ',pind,') => adjusted to lowest level, (pind=1.)'
 1079 
                    write(*,'(a,f7.3,a)') ' WARNING: pressure too high (pind = ',pind,') => adjusted to lowest level, (pind=1.)'
1076 
                    print*,'(x0,y0,p0)=',x0,y0,p0
 1080 
                    print*,'(x0,y0,p0)=',x0,y0,p0
1077 
                    pind = 1.
 1081 
                    pind = 1.
1078 
                 elseif ( err_c2.eq.10 ) then
 1082 
                 elseif ( err_c2.eq.10 ) then
1079 
                    print*,' WARNING: more pressures too high -> adjusted to lowest level '
 1083 
                    print*,' WARNING: more pressures too high -> adjusted to lowest level '
1080 
                    pind = 1.
 1084 
                    pind = 1.
Line 1429... Line 1433...
1429 
!            if ( (xind.ge.1.).and.(xind.le.real(nx)).and. &
 1433 
!            if ( (xind.ge.1.).and.(xind.le.real(nx)).and. &
1430 
!                 (yind.ge.1.).and.(yind.le.real(ny)).and. &
 1434 
!                 (yind.ge.1.).and.(yind.le.real(ny)).and. &
1431 
!                 (pind.ge.1.).and.(pind.le.real(nz)) ) then
 1435 
!                 (pind.ge.1.).and.(pind.le.real(nz)) ) then
1432 
!
 1436 
!
1433 
            ! Do the interpolation: everthing is ok
 1437 
            ! Do the interpolation: everthing is ok
- 
 
 1438 
            if ( ( shift_dir(i).ne.'PMIN' ).and.( shift_dir(i).ne.'PMAX' ) ) then
1434 
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
 1439 
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
1435 
 
 1440
- 
 
 1441 
            elseif  ( shift_dir(i).eq.'PMIN' ) then
- 
 
 1442 
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1443
- 
 
 1444 
               ! Handle > and < operator
- 
 
 1445 
               if ( (shift_rel(i).eq.'>').and.(f0.gt.shift_val(i)) ) then
- 
 
 1446 
                 do while ( (f0.gt.shift_val(i)).and.(pind.lt.nz) )
- 
 
 1447 
                     pind = pind + 0.1
- 
 
 1448 
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1449 
                 enddo
- 
 
 1450 
                 f0 = int_index4(p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1451
- 
 
 1452 
               elseif ( (shift_rel(i).eq.'<').and.(f0.lt.shift_val(i)) ) then
- 
 
 1453 
                 do while ( (f0.lt.shift_val(i)).and.(pind.lt.nz) )
- 
 
 1454 
                     pind = pind + 0.1
- 
 
 1455 
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1456 
                 enddo
- 
 
 1457 
                 f0 = int_index4(p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1458
- 
 
 1459 
              else
- 
 
 1460 
                 f0 = mdv
- 
 
 1461 
              endif
- 
 
 1462
- 
 
 1463
- 
 
 1464 
           elseif  ( shift_dir(i).eq.'PMAX' ) then
- 
 
 1465 
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1466
- 
 
 1467 
               ! Handle > and < operator
- 
 
 1468 
               if ( (shift_rel(i).eq.'>').and.(f0.gt.shift_val(i)) ) then
- 
 
 1469 
                 do while ( (f0.gt.shift_val(i)).and.(pind.gt.1) )
- 
 
 1470 
                     pind = pind - 0.1
- 
 
 1471 
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1472 
                 enddo
- 
 
 1473 
                 f0 = int_index4(p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1474
- 
 
 1475 
               elseif ( (shift_rel(i).eq.'<').and.(f0.lt.shift_val(i)) ) then
- 
 
 1476 
                 do while ( (f0.lt.shift_val(i)).and.(pind.gt.1) )
- 
 
 1477 
                     pind = pind - 0.1
- 
 
 1478 
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1479 
                 enddo
- 
 
 1480 
                 f0 = int_index4(p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1481
- 
 
 1482 
              else
- 
 
 1483 
                 f0 = mdv
- 
 
 1484 
              endif
- 
 
 1485
- 
 
 1486 
          endif
- 
 
 1487 
 
1436 
!            ! Check if pressure is outside, but rest okay: adjust to lowest or highest level
 1488 
!            ! Check if pressure is outside, but rest okay: adjust to lowest or highest level
1437 
!            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
 1489 
!            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
1438 
!               if ( pind.gt.nz ) then ! pressure too low, index too high
 1490 
!               if ( pind.gt.nz ) then ! pressure too low, index too high
1439 
!                 pind = real(nz)
 1491 
!                 pind = real(nz)
1440 
!                 print*,' Warning: pressure too low -> adjusted to highest level, pind=nz.'
 1492 
!                 print*,' Warning: pressure too low -> adjusted to highest level, pind=nz.'
Line 1464... Line 1516...
1464 
!            endif
 1516 
!            endif
1465 
 
 1517 
 
1466 
            ! ------------------------ end NORMAL mode -------------------------------
 1518 
            ! ------------------------ end NORMAL mode -------------------------------
1467 
            endif ! end if nearest case
 1519 
            endif ! end if nearest case
1468 
 
 1520 
 
1469 
           ! Exit for loop over all trajectories and times -save interpolated value
 1521 
           ! Exit for loop over all times -save interpolated value
1470 
 109        continue
 1522 
 109        continue
1471 
 
 1523 
 
1472 
            ! Save the new field
 1524 
            ! Save the new field
1473 
            if ( abs(f0-mdv).gt.eps) then
 1525 
            if ( abs(f0-mdv).gt.eps) then
1474 
               traint(k,j,ncol+i) = f0
 1526 
               traint(k,j,ncol+i) = f0
Line 1976... Line 2028...
1976 
 
 2028 
 
1977 
  ! Set the variable names for output trajectory
 2029 
  ! Set the variable names for output trajectory
1978 
  do i=1,ncol+ntrace0
 2030 
  do i=1,ncol+ntrace0
1979 
     varsout(i)      = varsint(i)
 2031 
        varsout(i)      = varsint(i)
1980 
  enddo
 2032 
  enddo
- 
 
 2033 
  do i=1,ntrace0
- 
 
 2034 
     if ( (shift_dir(i).eq.'PMIN').or.(shift_dir(i).eq.'PMAX') ) then
- 
 
 2035 
        varsout(ncol+i) = trim(tvar(i))//':'//trim(shift_dir(i))
- 
 
 2036 
     endif
- 
 
 2037 
  enddo
1981 
 
 2038
1982 
  ! Write trajectories
 2039 
  ! Write trajectories
1983 
  call wopen_tra(fod,outfile,ntra,ntim,ncol+ntrace0, &
 2040 
  call wopen_tra(fod,outfile,ntra,ntim,ncol+ntrace0, &
1984 
       reftime,varsout,outmode)
 2041 
       reftime,varsout,outmode)
1985 
  call write_tra(fod,traout ,ntra,ntim,ncol+ntrace0,outmode)
 2042 
  call write_tra(fod,traout ,ntra,ntim,ncol+ntrace0,outmode)
Line 2093... Line 2150...
2093 
 
 2150 
 
2094 
! ------------------------------------------------------------------
 2151 
! ------------------------------------------------------------------
2095 
! Split the variable name into name, shift and direction
 2152 
! Split the variable name into name, shift and direction
2096 
! ------------------------------------------------------------------
 2153 
! ------------------------------------------------------------------
2097 
 
 2154 
 
2098 
subroutine splitvar (tvar,shift_val,shift_dir)
 2155 
subroutine splitvar (tvar,shift_val,shift_dir,shift_rel)
2099 
 
 2156 
 
2100 
  implicit none
 2157 
  implicit none
2101 
 
 2158 
 
2102 
  ! Declaration of subroutine parameters
 2159 
  ! Declaration of subroutine parameters
2103 
  character(len=80) :: tvar
 2160 
  character(len=80) :: tvar
2104 
  real :: shift_val
 2161 
  real :: shift_val
2105 
  character(len=80) :: shift_dir
 2162 
  character(len=80) :: shift_dir
- 
 
 2163 
  character(len=80) :: shift_rel
2106 
 
 2164 
 
2107 
  ! Auxiliary variables
 2165 
  ! Auxiliary variables
2108 
  integer :: i,j
 2166 
  integer :: i,j
2109 
  integer :: icolon,inumber
 2167 
  integer :: icolon,inumber,irelator
2110 
  character(len=80) :: name
 2168 
  character(len=80) :: name
2111 
  character :: ch
 2169 
  character :: ch
2112 
  integer      isabsval
 2170 
  integer      isabsval
2113 
 
 2171 
 
2114 
  ! Save variable name
 2172 
  ! Save variable name
2115 
  name = tvar
 2173 
  name = tvar
- 
 
 2174 
  shift_rel = 'nil'
- 
 
 2175 
  shift_dir = 'nil'
2116 
 
 2176 
 
2117 
  ! Search for colon
 2177 
  ! Search for colon
2118 
  icolon=0
 2178 
  icolon=0
2119 
  do i=1,80
 2179 
  do i=1,80
2120 
     if ( (name(i:i).eq.':').and.(icolon.ne.0) ) goto 100
 2180 
     if ( (name(i:i).eq.':').and.(icolon.ne.0) ) goto 100
Line 2137... Line 2197...
2137 
           inumber = i
 2197 
           inumber = i
2138 
           exit
 2198 
           exit
2139 
        endif
 2199 
        endif
2140 
     enddo
 2200 
     enddo
2141 
 
 2201 
 
- 
 
 2202 
     ! If the variable name is e.g. PMIN:UMF>0, the variable to be read
- 
 
 2203 
     ! is UMF, the value is 0, and the direction is 'PMIN'
- 
 
 2204 
     shift_rel = 'nil'
- 
 
 2205 
     if ( (tvar.eq.'PMIN').or.(tvar.eq.'PMAX') ) then
- 
 
 2206 
         shift_dir = tvar
- 
 
 2207 
         irelator = 0
- 
 
 2208 
         do i=icolon+1,80
- 
 
 2209 
            ch = name(i:i)
- 
 
 2210 
            if ( (ch.eq.'>').or.(ch.eq.'<') ) then
- 
 
 2211 
               irelator = i
- 
 
 2212 
            endif
- 
 
 2213 
         enddo
- 
 
 2214 
         if ( irelator.eq.0 ) then
- 
 
 2215 
             print*,' ERROR: dont know how to interpret ',trim(name)
- 
 
 2216 
             stop
- 
 
 2217 
         endif
- 
 
 2218 
         tvar = name(icolon+1:irelator-1)
- 
 
 2219 
         read(name( (irelator+1):80 ),*) shift_val
- 
 
 2220 
         shift_rel = name(irelator:irelator)
- 
 
 2221
- 
 
 2222 
         goto 90
- 
 
 2223
- 
 
 2224 
     endif
- 
 
 2225
2142 
     ! Get the number
 2226 
     ! Get the number
2143 
     read(name( (icolon+1):(inumber-1) ),*) shift_val
 2227 
     read(name( (icolon+1):(inumber-1) ),*) shift_val
2144 
 
 2228 
 
2145 
     ! Decide whether it is a shift relatiev to trajectory or absolute value
 2229 
     ! Decide whether it is a shift relatiev to trajectory or absolute value
2146 
     ! If the number starts with + or -, it is relative to the trajectory
 2230 
     ! If the number starts with + or -, it is relative to the trajectory
Line 2161... Line 2245...
2161 
     shift_dir = 'nil'
 2245 
     shift_dir = 'nil'
2162 
     shift_val = 0.
 2246 
     shift_val = 0.
2163 
 
 2247 
 
2164 
  endif
 2248 
  endif
2165 
 
 2249 
 
- 
 
 2250 
 90 continue
- 
 
 2251
2166 
  return
 2252 
    return
2167 
 
 2253 
 
2168 
  ! Error handling
 2254 
  ! Error handling
2169 
 100   continue
 2255 
 100   continue
2170 
 
 2256