Subversion Repositories lagranto.ecmwf

Rev

Rev 13 | Rev 27 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
      PROGRAM caltra
2
 
3
C     ********************************************************************
4
C     *                                                                  *
5
C     * Calculates trajectories                                          *
6
C     *                                                                  *
7
C     *	Heini Wernli	   first version:	April 1993               *
8
C     * Michael Sprenger   major upgrade:       2008-2009                *
9
C     *                                                                  *
10
C     ********************************************************************
11
 
12
      implicit none      
13
 
14
c     --------------------------------------------------------------------
15
c     Declaration of parameters
16
c     --------------------------------------------------------------------
17
 
18
c     Maximum number of levels for input files
19
      integer   nlevmax
20
      parameter	(nlevmax=100)
21
 
22
c     Maximum number of input files (dates, length of trajectories)
23
      integer   ndatmax
24
      parameter	(ndatmax=500)
25
 
26
c     Numerical epsilon (for float comparison)
27
      real      eps
28
      parameter (eps=0.001)
29
 
30
c     Distance in m between 2 lat circles 
31
      real	deltay
32
      parameter	(deltay=1.112E5)
33
 
34
c     Numerical method for the integration (0=iterative Euler, 1=Runge-Kutta)
35
      integer   imethod
36
      parameter (imethod=1)
37
 
38
c     Number of iterations for iterative Euler scheme
39
      integer   numit
40
      parameter (numit=3)
41
 
42
c     Input and output format for trajectories (see iotra.f)
43
      integer   inpmode
44
      integer   outmode
45
 
19 michaesp 46
c     Filename prefix for primary and secondary files (typically 'P' and 'S')
3 michaesp 47
      character*1 prefix
48
      parameter   (prefix='P')
19 michaesp 49
      character*1 srefix
50
      parameter   (srefix='S')
3 michaesp 51
 
19 michaesp 52
c     Physical constants - needed to compute potential temperature
53
      real      rdcp,tzero
54
      data      rdcp,tzero /0.286,273.15/
55
 
3 michaesp 56
c     --------------------------------------------------------------------
57
c     Declaration of variables
58
c     --------------------------------------------------------------------
59
 
60
c     Input parameters
61
      integer                                fbflag          ! Flag for forward/backward mode
62
      integer                                numdat          ! Number of input files
63
      character*11                           dat(ndatmax)    ! Dates of input files
64
      real                                   timeinc         ! Time increment between input files
65
      real                                   per             ! Periodicity (=0 if none)
66
      integer                                ntra            ! Number of trajectories
67
      character*80                           cdfname         ! Name of output files
68
      real                                   ts              ! Time step
69
      real                                   tst,ten         ! Shift of start and end time relative to first data file
70
      integer                                deltout         ! Output time interval (in minutes)
71
      integer                                jflag           ! Jump flag (if =1 ground-touching trajectories reenter atmosphere)
72
      real                                   wfactor         ! Factor for vertical velocity field
73
      character*80                           strname         ! File with start positions
74
      character*80                           timecheck       ! Either 'yes' or 'no'
19 michaesp 75
      character*80                           isen            ! Isentropic trajectories ('yes' or 'no')
76
      integer                                thons           ! Isentropic mode: is TH availanle on S
77
      character*80                           modlev          ! 2D (model level) trajectories ('yes' or 'no')
3 michaesp 78
 
79
c     Trajectories
80
      integer                                ncol            ! Number of columns for insput trajectories
81
      real,allocatable, dimension (:,:,:) :: trainp          ! Input start coordinates (ntra,1,ncol)
82
      real,allocatable, dimension (:,:,:) :: traout          ! Output trajectories (ntra,ntim,4)
83
      integer                                reftime(6)      ! Reference date
84
      character*80                           vars(200)       ! Field names
85
      real,allocatable, dimension (:)     :: xx0,yy0,pp0     ! Position of air parcels
86
      integer,allocatable, dimension (:)  :: leftflag        ! Flag for domain-leaving
87
      real                                   xx1,yy1,pp1     ! Updated position of air parcel
88
      integer                                leftcount       ! Number of domain leaving trajectories
89
      integer                                ntim            ! Number of output time steps
19 michaesp 90
      real,allocatable, dimension (:)     :: theta           ! Potential temperature for isentropic trajectories
91
      real,allocatable, dimension (:)     :: zindex          ! Vertical index for model-level (2D) trajectories
3 michaesp 92
 
93
c     Meteorological fields
94
      real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
95
      real,allocatable, dimension (:)     :: uut0,uut1       ! Zonal wind
96
      real,allocatable, dimension (:)     :: vvt0,vvt1       ! Meridional wind
97
      real,allocatable, dimension (:)     :: wwt0,wwt1       ! Vertical wind
98
      real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure 
19 michaesp 99
      real,allocatable, dimension (:)     :: tht0,tht1       ! 3d potential temperature
100
      real,allocatable, dimension (:)     :: sth0,sth1       ! Surface potential temperature
3 michaesp 101
 
102
c     Grid description
103
      real                                   pollon,pollat   ! Longitude/latitude of pole
104
      real                                   ak(nlevmax)     ! Vertical layers and levels
105
      real                                   bk(nlevmax) 
106
      real                                   xmin,xmax       ! Zonal grid extension
107
      real                                   ymin,ymax       ! Meridional grid extension
108
      integer                                nx,ny,nz        ! Grid dimensions
109
      real                                   dx,dy           ! Horizontal grid resolution
110
      integer                                hem             ! Flag for hemispheric domain
111
      real                                   mdv             ! Missing data value
112
 
113
c     Auxiliary variables                 
114
      real                                   delta,rd
115
      integer	                             itm,iloop,i,j,k,filo,lalo
116
      integer                                ierr,stat
117
      integer                                cdfid,fid
9 michaesp 118
      real	                                 tstart,time0,time1,time
3 michaesp 119
      real                                   reltpos0,reltpos1
120
      real                                   xind,yind,pind,pp,sp,stagz
121
      character*80                           filename,varname
122
      integer                                reftmp(6)
123
      character                              ch
124
      real                                   frac,tload
125
      integer                                itim
126
      integer                                wstep
19 michaesp 127
      real                                   x1,y1,p1
128
      real                                   thetamin,thetamax
129
      real                                   zindexmin,zindexmax
3 michaesp 130
 
9 michaesp 131
c     Externals
132
      real                                   int_index4
133
      external                               int_index4
134
 
3 michaesp 135
c     --------------------------------------------------------------------
136
c     Start of program, Read parameters
137
c     --------------------------------------------------------------------
138
 
139
c     Write start message
140
      print*,'========================================================='
141
      print*,'              *** START OF PROGRAM CALTRA ***'
142
      print*
143
 
144
c     Open the parameter file
145
      open(9,file='caltra.param')
146
 
147
c     Read flag for forward/backward mode (fbflag)
148
      read(9,*) fbflag
149
 
150
c     Read number of input files (numdat)
151
      read(9,*) numdat
152
      if (numdat.gt.ndatmax) then
153
        print*,' ERROR: too many input files ',numdat,ndatmax
154
        goto 993
155
      endif
156
 
157
c     Read list of input dates (dat, sort depending on forward/backward mode)
158
      if (fbflag.eq.1) then
159
        do itm=1,numdat
160
          read(9,'(a11)') dat(itm)
161
        enddo
162
      else
163
        do itm=numdat,1,-1
164
          read(9,'(a11)') dat(itm)
165
        enddo
166
      endif
167
 
168
c     Read time increment between input files (timeinc)
169
      read(9,*) timeinc
170
 
171
C     Read if data domain is periodic and its periodicity
172
      read(9,*) per
173
 
174
c     Read the number of trajectories and name of position file
175
      read(9,*) strname
176
      read(9,*) ntra
177
      read(9,*) ncol 
178
      if (ntra.eq.0) goto 991
179
 
180
C     Read the name of the output trajectory file and set the constants file
181
      read(9,*) cdfname
182
 
183
C     Read the timestep for trajectory calculation (convert from minutes to hours)
184
      read(9,*) ts
185
      ts=ts/60.        
186
 
187
C     Read shift of start and end time relative to first data file
188
      read(9,*) tst
189
      read(9,*) ten
190
 
191
C     Read output time interval (in minutes)
192
      read(9,*) deltout
193
 
194
C     Read jumpflag (if =1 ground-touching trajectories reenter the atmosphere)
195
      read(9,*) jflag
196
 
197
C     Read factor for vertical velocity field
198
      read(9,*) wfactor
199
 
200
c     Read the reference time and the time range
201
      read(9,*) reftime(1)                  ! year
202
      read(9,*) reftime(2)                  ! month
203
      read(9,*) reftime(3)                  ! day
204
      read(9,*) reftime(4)                  ! hour
205
      read(9,*) reftime(5)                  ! min
206
      read(9,*) reftime(6)                  ! time range (in min)
207
 
208
c     Read flag for 'no time check'
209
      read(9,*) timecheck
210
 
19 michaesp 211
c     Read flag for isentropic trajectories
212
      read(9,*) isen, thons
213
 
214
c     Read flag for model-level trajectories (2D mode)
215
      read(9,*) modlev
216
 
3 michaesp 217
c     Close the input file
218
      close(9)
219
 
220
c     Calculate the number of output time steps
221
      ntim = abs(reftime(6)/deltout) + 1
222
 
223
c     Set the formats of the input and output files
224
      call mode_tra(inpmode,strname)
225
      call mode_tra(outmode,cdfname)
226
      if (outmode.eq.-1) outmode=1
227
 
228
c     Write some status information
229
      print*,'---- INPUT PARAMETERS -----------------------------------'
230
      print* 
231
      print*,'  Forward/Backward       : ',fbflag
232
      print*,'  #input files           : ',numdat
233
      print*,'  First/last input file  : ',trim(dat(1)),' ... ',
234
     >                                     trim(dat(numdat))
235
      print*,'  time increment         : ',timeinc
236
      print*,'  Output file            : ',trim(cdfname)
237
      print*,'  Time step (min)        : ',60.*ts
238
      write(*,'(a27,f7.2,f7.2)') '   Time shift (start,end) : ',tst,ten
239
      print*,'  Output time interval   : ',deltout
240
      print*,'  Jump flag              : ',jflag
241
      print*,'  Vertical wind (scale)  : ',wfactor
242
      print*,'  Trajectory pos file    : ',trim(strname)
243
      print*,'  # of trajectories      : ',ntra
244
      print*,'  # of output timesteps  : ',ntim
245
      if ( inpmode.eq.-1) then
246
         print*,'  Input format           : (lon,lat,p)-list'
247
      else
248
         print*,'  Input format           : ',inpmode
249
      endif
250
      print*,'  Output format          : ',outmode
251
      print*,'  Periodicity            : ',per
252
      print*,'  Time check             : ',trim(timecheck)
19 michaesp 253
      print*,'  Isentropic trajectories: ',trim(isen),thons
254
      print*,'  Model-level trajs (2D) : ',trim(modlev)
3 michaesp 255
      print*
256
 
19 michaesp 257
      if ( (isen.eq.'yes').and.(modlev.eq.'yes') ) then
258
         print*,
259
     >    '  WARNING: isentropic and 2D mode chosen -> 2D accepted'
260
         print*
261
         isen = 'no'
262
      endif
263
 
264
c     Init missing data value
265
      mdv = -999.
266
 
3 michaesp 267
      print*,'---- FIXED NUMERICAL PARAMETERS -------------------------'
268
      print*
269
      print*,'  Numerical scheme       : ',imethod
270
      print*,'  Number of iterations   : ',numit
271
      print*,'  Filename prefix        : ',prefix
272
      print*,'  Missing data value     : ',mdv
273
      print*
274
 
275
c     --------------------------------------------------------------------
276
c     Read grid parameters, checks and allocate memory
277
c     --------------------------------------------------------------------
278
 
279
c     Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)
280
c     The negative <-fid> of the file identifier is used as a flag for parameter retrieval  
281
      filename = prefix//dat(1)
282
      varname  = 'U'
283
      nx       = 1
284
      ny       = 1
285
      nz       = 1
286
      tload    = -tst
287
      call input_open (fid,filename)
288
      call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
289
     >                 tload,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)
290
      call input_close(fid)
291
 
292
C     Check if the number of levels is too large
293
      if (nz.gt.nlevmax) goto 993
294
 
295
C     Set logical flag for periodic data set (hemispheric or not)
9 michaesp 296
      hem = 0
3 michaesp 297
      if (per.eq.0.) then
298
         delta=xmax-xmin-360.
299
         if (abs(delta+dx).lt.eps) then               ! Program aborts: arrays must be closed
300
            goto 992
301
        else if (abs(delta).lt.eps) then              ! Periodic and hemispheric
302
           hem=1
303
           per=360.
304
        endif
305
      else                                            ! Periodic and hemispheric
306
         hem=1
307
      endif
308
 
309
C     Allocate memory for some meteorological arrays
310
      allocate(spt0(nx*ny),stat=stat)
311
      if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface pressure
312
      allocate(spt1(nx*ny),stat=stat)
313
      if (stat.ne.0) print*,'*** error allocating array spt1 ***'
314
      allocate(uut0(nx*ny*nz),stat=stat)
315
      if (stat.ne.0) print*,'*** error allocating array uut0 ***'   ! Zonal wind
316
      allocate(uut1(nx*ny*nz),stat=stat)
317
      if (stat.ne.0) print*,'*** error allocating array uut1 ***'
318
      allocate(vvt0(nx*ny*nz),stat=stat)
319
      if (stat.ne.0) print*,'*** error allocating array vvt0 ***'   ! Meridional wind
320
      allocate(vvt1(nx*ny*nz),stat=stat)
321
      if (stat.ne.0) print*,'*** error allocating array vvt1 ***'
322
      allocate(wwt0(nx*ny*nz),stat=stat)
323
      if (stat.ne.0) print*,'*** error allocating array wwt0 ***'   ! Vertical wind
324
      allocate(wwt1(nx*ny*nz),stat=stat)
325
      if (stat.ne.0) print*,'*** error allocating array wwt1 ***'
326
      allocate(p3t0(nx*ny*nz),stat=stat)
327
      if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Pressure
328
      allocate(p3t1(nx*ny*nz),stat=stat)
329
      if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
19 michaesp 330
      allocate(tht0(nx*ny*nz),stat=stat)
331
      if (stat.ne.0) print*,'*** error allocating array tht0 ***'   ! Potential temperature
332
      allocate(tht1(nx*ny*nz),stat=stat)
333
      if (stat.ne.0) print*,'*** error allocating array tht1 ***'
334
      allocate(sth0(nx*ny),stat=stat)
335
      if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface potential temperature
336
      allocate(sth1(nx*ny),stat=stat)
337
      if (stat.ne.0) print*,'*** error allocating array sth1 ***'
3 michaesp 338
 
339
C     Get memory for trajectory arrays
340
      allocate(trainp(ntra,1,ncol),stat=stat)
341
      if (stat.ne.0) print*,'*** error allocating array trainp   ***' ! Input start coordinates
342
      allocate(traout(ntra,ntim,4),stat=stat)
343
      if (stat.ne.0) print*,'*** error allocating array traout   ***' ! Output trajectories
344
      allocate(xx0(ntra),stat=stat)
345
      if (stat.ne.0) print*,'*** error allocating array xx0      ***' ! X position (longitude)
346
      allocate(yy0(ntra),stat=stat)
347
      if (stat.ne.0) print*,'*** error allocating array yy0      ***' ! Y position (latitude)
348
      allocate(pp0(ntra),stat=stat)
349
      if (stat.ne.0) print*,'*** error allocating array pp0      ***' ! Pressure
350
      allocate(leftflag(ntra),stat=stat)
351
      if (stat.ne.0) print*,'*** error allocating array leftflag ***' ! Leaving-domain flag
19 michaesp 352
      allocate(theta(ntra),stat=stat)
353
      if (stat.ne.0) print*,'*** error allocating array theta ***'    ! Potential temperature for isentropic trajectories
354
      allocate(zindex(ntra),stat=stat)
355
      if (stat.ne.0) print*,'*** error allocating array kind ***'     ! Vertical index for model-level trajectories
3 michaesp 356
 
357
c     Write some status information
358
      print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
359
      print*
360
      print*,'  xmin,xmax     : ',xmin,xmax
361
      print*,'  ymin,ymax     : ',ymin,ymax
362
      print*,'  dx,dy         : ',dx,dy
363
      print*,'  pollon,pollat : ',pollon,pollat
364
      print*,'  nx,ny,nz      : ',nx,ny,nz
365
      print*,'  per, hem      : ',per,hem
366
      print*
367
 
368
c     --------------------------------------------------------------------
369
c     Initialize the trajectory calculation
370
c     --------------------------------------------------------------------
371
 
372
c     Read start coordinates from file - Format (lon,lat,lev)
373
      if (inpmode.eq.-1) then
374
         open(fid,file=strname)
375
          do i=1,ntra
376
             read(fid,*) xx0(i),yy0(i),pp0(i)
377
          enddo
378
         close(fid)
379
 
380
c     Read start coordinates from trajectory file - check consistency of ref time
381
      else
382
         call ropen_tra(cdfid,strname,ntra,1,ncol,reftmp,vars,inpmode)
383
         call read_tra (cdfid,trainp,ntra,1,ncol,inpmode)
384
         do i=1,ntra
385
            time   = trainp(i,1,1)
386
            xx0(i) = trainp(i,1,2) 
387
            yy0(i) = trainp(i,1,3) 
388
            pp0(i) = trainp(i,1,4) 
389
         enddo
390
         call close_tra(cdfid,inpmode)
391
 
392
         if ( ( reftime(1).ne.reftmp(1) ).or.
393
     >        ( reftime(2).ne.reftmp(2) ).or.
394
     >        ( reftime(3).ne.reftmp(3) ).or.
395
     >        ( reftime(4).ne.reftmp(4) ).or.
396
     >        ( reftime(5).ne.reftmp(5) ) )
397
     >   then
398
            print*,' WARNING: Inconsistent reference times'
399
            write(*,'(5i8)') (reftime(i),i=1,5)
400
            write(*,'(5i8)') (reftmp (i),i=1,5)
401
            print*,'Enter a key to proceed...'
402
            stop
403
         endif
404
      endif
405
 
406
c     Set sign of time range
407
      reftime(6) = fbflag * reftime(6)
408
 
409
c     Write some status information
410
      print*,'---- REFERENCE DATE---------- ---------------------------'
411
      print*
412
      print*,' Reference time (year)  :',reftime(1)
413
      print*,'                (month) :',reftime(2)
414
      print*,'                (day)   :',reftime(3)
415
      print*,'                (hour)  :',reftime(4)
416
      print*,'                (min)   :',reftime(5)
417
      print*,' Time range             :',reftime(6),' min'
418
      print*
419
 
420
C     Save starting positions 
421
      itim = 1
422
      do i=1,ntra
423
         traout(i,itim,1) = 0.
424
         traout(i,itim,2) = xx0(i)
425
         traout(i,itim,3) = yy0(i)
426
         traout(i,itim,4) = pp0(i)
427
      enddo
9 michaesp 428
 
3 michaesp 429
c     Init the flag and the counter for trajectories leaving the domain
430
      leftcount=0
431
      do i=1,ntra
432
         leftflag(i)=0
433
      enddo
434
 
435
C     Convert time shifts <tst,ten> from <hh.mm> into fractional time
436
      call hhmm2frac(tst,frac)
437
      tst = frac
438
      call hhmm2frac(ten,frac)
439
      ten = frac
440
 
19 michaesp 441
c     Get 3D and surface pressure from first data file 
9 michaesp 442
      filename = prefix//dat(1)
443
      call input_open (fid,filename)
19 michaesp 444
      if ( modlev.eq.'no' ) then 
445
         varname = 'P'
446
         call input_grid
447
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
448
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
449
      else  
450
         varname = 'P.ML'
451
         call input_grid
452
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
453
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
454
      endif
9 michaesp 455
      call input_close(fid)
19 michaesp 456
 
457
c     Check that all starting positions are above topography    
9 michaesp 458
      do i=1,ntra
459
 
460
C       Interpolate surface pressure to actual position (from first input file)
461
        x1 = xx0(i)
462
        y1 = yy0(i)
463
        call get_index4 (xind,yind,pind,x1,y1,1050.,0.,
464
     >                   p3t1,p3t1,spt1,spt1,3,
465
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
466
        sp = int_index4 (spt1,spt1,nx,ny,1,xind,yind,1.,0.,mdv)
467
 
468
c       Decide whether to keep the trajectory
469
        if ( pp0(i).gt.sp ) then
470
            write(*,'(a30,3f10.2)')
471
     >               'WARNING: starting point below topography ',
472
     >               xx0(i),yy0(i),pp0(i)
473
            leftflag(i) = 1
474
        endif
475
 
476
      enddo
477
 
19 michaesp 478
c     Special handling for isentropic trajectories - read potential
479
c     temperature from S file or calculate it based on temperature and
480
c     pressure from P file; then, calculate for each trajectory its
481
c     potential temperature - which will stay fixed over time
482
      if ( isen.eq.'yes' ) then
9 michaesp 483
 
19 michaesp 484
c         Get potential temperature from S file
485
          if ( thons.eq.1 ) then
486
              filename = srefix//dat(1)
487
              print*,' TH <- ',trim(filename)
488
              call input_open (fid,filename)
489
              varname='TH'
490
              call input_wind
491
     >            (fid,varname,tht1,tload,stagz,mdv,
492
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
493
              call input_close(fid)
494
 
495
c         Calculate potential temperature from P file
496
          else
497
              filename = prefix//dat(1)
498
              print*,' TH = T * (1000/P)^RDCP <- ',trim(filename)
499
              call input_open (fid,filename)
500
              varname='T'
501
              call input_wind
502
     >            (fid,varname,tht1,tload,stagz,mdv,
503
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
504
              call input_close(fid)
505
              do i=1,nx*ny*nz
506
                if (tht1(i).lt.100.) then
507
                   tht1(i)=(tht1(i)+tzero)*( (1000./p3t1(i))**rdcp )
508
                else
509
                   tht1(i)=tht1(i)*( (1000./p3t1(i))**rdcp )
510
                endif
511
              enddo
512
          endif
513
 
514
c         Take surface potential temperature from lowest level
515
          do i=1,nx*ny
516
            sth1(i) = tht1(i)
517
          enddo
518
 
519
c         Calculate now the potential temperature of all trajectories
520
          do i=1,ntra
521
 
522
             x1 = xx0(i)
523
             y1 = yy0(i)
524
             p1 = pp0(i)
525
             call get_index4 (xind,yind,pind,x1,y1,p1,0.,
526
     >                 p3t1,p3t1,spt1,spt1,3,
527
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
528
             theta(i) =
529
     >        int_index4(tht1,tht1,nx,ny,nz,xind,yind,pind,0.,mdv)
530
 
531
          enddo
532
 
533
c         Write info about theta range of starting positions
534
          thetamin = theta(1)
535
          thetamax = theta(1)
536
          do i=2,ntra
537
            if ( theta(i).gt.thetamax ) thetamax = theta(i)
538
            if ( theta(i).lt.thetamin ) thetamin = theta(i)
539
          enddo
540
 
541
c         Write some status information
542
          print*
543
          print*,
544
     >     '---- THETA RANGE OF ISENTROPIC TRAJECTORIES -------------'
545
          print*
546
          print*,' Theta(min)             :',thetamin
547
          print*,' Theta(max)             :',thetamax
548
 
549
      endif
550
 
551
c     Special handling for model-level (2D) trajectories - get the
552
c     vertical index for each trajectory - which will remain fixed
553
      if ( modlev.eq.'yes' ) then
554
          do i=1,ntra
555
             x1 = xx0(i)
556
             y1 = yy0(i)
557
             p1 = pp0(i)
558
             call get_index4 (xind,yind,pind,x1,y1,p1,0.,
559
     >                 p3t1,p3t1,spt1,spt1,3,
560
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
561
             zindex(i) = pind
562
          enddo
563
 
564
          do i=1,nz
565
            print*,i,p3t1(189+(141-1)*nx+(i-1)*nx*ny)
566
          enddo
567
             print*,x1,y1,p1
568
             print*,xind,yind,pind
569
 
570
c         Write info about zindex range of starting positions
571
          zindexmin = zindex(1)
572
          zindexmax = zindex(1)
573
          do i=2,ntra
574
            if ( zindex(i).gt.zindexmax ) zindexmax = zindex(i)
575
            if ( zindex(i).lt.zindexmin ) zindexmin = zindex(i)
576
          enddo
577
 
578
c         Write some status information
579
          print*
580
          print*,
581
     >     '---- INDEX RANGE OF MODEL-LEVEL TRAJECTORIES ------------'
582
          print*
583
          print*,' Zindex(min)            :',zindexmin
584
          print*,' Zindex(max)            :',zindexmax
585
 
586
 
587
      endif
588
 
3 michaesp 589
c     -----------------------------------------------------------------------
590
c     Loop to calculate trajectories
591
c     -----------------------------------------------------------------------
592
 
593
c     Write some status information
9 michaesp 594
      print*
3 michaesp 595
      print*,'---- TRAJECTORIES ----------- ---------------------------'
596
      print*    
597
 
598
C     Set the time for the first data file (depending on forward/backward mode)
599
      if (fbflag.eq.1) then
600
        tstart = -tst
601
      else
602
        tstart = tst
603
      endif
604
 
605
c     Set the minute counter for output
606
      wstep = 0
607
 
608
c     Read wind fields and vertical grid from first file
609
      filename = prefix//dat(1)
610
 
611
      call frac2hhmm(tstart,tload)
612
 
19 michaesp 613
      write(*,'(a16,a20,f9.2)') '  (file,time) : ',
3 michaesp 614
     >                       trim(filename),tload
615
 
616
      call input_open (fid,filename)
19 michaesp 617
 
3 michaesp 618
      varname='U'                                      ! U
619
      call input_wind 
620
     >    (fid,varname,uut1,tload,stagz,mdv,
621
     >     xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
19 michaesp 622
 
3 michaesp 623
      varname='V'                                      ! V
624
      call input_wind 
625
     >    (fid,varname,vvt1,tload,stagz,mdv,
626
     >     xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
19 michaesp 627
 
628
      if ( (modlev.eq.'no').and.(isen.eq.'no') ) then
629
         varname='OMEGA'                                  ! OMEGA
630
         call input_wind
631
     >       (fid,varname,wwt1,tload,stagz,mdv,
632
     >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
633
      endif
634
 
635
      if ( modlev.eq.'no' ) then
636
         call input_grid                                  ! GRID - AK,BK -> P
637
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
638
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
639
      else
640
         varname='P.ML'                                   ! GRID - P,PS
641
         call input_grid                                  !
642
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
643
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
644
      endif
645
 
3 michaesp 646
      call input_close(fid)
19 michaesp 647
 
648
c     Special handling for isentropic trajectories - read potential
649
c     temperature from S file or calculate it based on temperature and
650
c     pressure from P file
651
      if ( isen.eq.'yes' ) then
652
 
653
c         Get potential temperature from S file
654
          if ( thons.eq.1 ) then
655
              filename = srefix//dat(1)
656
              print*,' TH <- ',trim(filename)
657
              call input_open (fid,filename)
658
              varname='TH'
659
              call input_wind
660
     >            (fid,varname,tht1,tload,stagz,mdv,
661
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
662
              call input_close(fid)
663
 
664
c         Calculate potential temperature from P file
665
          else
666
              filename = prefix//dat(1)
667
              print*,' TH = T * (1000/P)^RDCP <- ',trim(filename)
668
              call input_open (fid,filename)
669
              varname='T'
670
              call input_wind
671
     >            (fid,varname,tht1,tload,stagz,mdv,
672
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
673
              call input_close(fid)
674
              do i=1,nx*ny*nz
675
                if (tht1(i).lt.100.) then
676
                   tht1(i)=(tht1(i)+tzero)*( (1000./p3t1(i))**rdcp )
677
                else
678
                   tht1(i)=tht1(i)*( (1000./p3t1(i))**rdcp )
679
                endif
680
              enddo
681
          endif
682
 
683
c         Take surface potential temperature from lowest level
684
          do i=1,nx*ny
685
            sth1(i) = tht1(i)
686
          enddo
687
      endif
688
 
3 michaesp 689
c     Loop over all input files (time step is <timeinc>)
690
      do itm=1,numdat-1
691
 
692
c       Calculate actual and next time
693
        time0 = tstart+real(itm-1)*timeinc*fbflag
694
        time1 = time0+timeinc*fbflag
695
 
696
c       Copy old velocities and pressure fields to new ones       
697
        do i=1,nx*ny*nz
698
           uut0(i)=uut1(i)
699
           vvt0(i)=vvt1(i)
700
           wwt0(i)=wwt1(i)
701
           p3t0(i)=p3t1(i)
19 michaesp 702
           tht0(i)=tht1(i)
3 michaesp 703
        enddo
704
        do i=1,nx*ny
705
           spt0(i)=spt1(i)
19 michaesp 706
           sth0(i)=sth1(i)
3 michaesp 707
        enddo
708
 
709
c       Read wind fields and surface pressure at next time
710
        filename = prefix//dat(itm+1)
711
 
712
        call frac2hhmm(time1,tload)
19 michaesp 713
        write(*,'(a16,a20,f9.2)') '  (file,time) : ',
3 michaesp 714
     >                          trim(filename),tload
715
 
716
        call input_open (fid,filename)
19 michaesp 717
 
3 michaesp 718
        varname='U'                                     ! U
719
        call input_wind 
720
     >       (fid,varname,uut1,tload,stagz,mdv,
721
     >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
19 michaesp 722
 
3 michaesp 723
        varname='V'                                     ! V
724
        call input_wind 
725
     >       (fid,varname,vvt1,tload,stagz,mdv,
726
     >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
19 michaesp 727
 
728
        if ( (modlev.eq.'no').and.(isen.eq.'no') ) then
729
           varname='OMEGA'                              ! OMEGA
730
           call input_wind
731
     >          (fid,varname,wwt1,tload,stagz,mdv,
732
     >           xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
733
        endif
734
 
735
        if ( modlev.eq.'no' ) then
736
           call input_grid                                  ! GRID - AK,NK -> P
737
     >          (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
738
     >           tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
739
        else
740
           varname='P.ML'                                   ! GRID - P,PS
741
           call input_grid                                  !
3 michaesp 742
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
743
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
19 michaesp 744
        endif
745
 
3 michaesp 746
        call input_close(fid)
19 michaesp 747
 
748
c     Special handling for isentropic trajectories - read potential
749
c     temperature from S file or calculate it based on temperature and
750
c     pressure from P file
751
      if ( isen.eq.'yes' ) then
752
 
753
c         Get TH from S file
754
          if ( thons.eq.1 ) then
755
              filename = srefix//dat(itm+1)
756
              print*,' TH <- ',trim(filename)
757
              call input_open (fid,filename)
758
              varname='TH'
759
              call input_wind
760
     >            (fid,varname,tht1,tload,stagz,mdv,
761
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
762
              call input_close(fid)
763
 
764
c         Calculate potential temperature from P file
765
          else
766
              filename = prefix//dat(itm+1)
767
              print*,' TH = T * (1000/P)^RDCP <- ',trim(filename)
768
              call input_open (fid,filename)
769
              varname='T'
770
              call input_wind
771
     >            (fid,varname,tht1,tload,stagz,mdv,
772
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
773
              call input_close(fid)
774
              do i=1,nx*ny*nz
775
                if (tht1(i).lt.100.) then
776
                   tht1(i)=(tht1(i)+tzero)*( (1000./p3t1(i))**rdcp )
777
                else
778
                   tht1(i)=tht1(i)*( (1000./p3t1(i))**rdcp )
779
                endif
780
              enddo
781
          endif
782
 
783
c         Take surface potential temperature from lowest level
784
          do i=1,nx*ny
785
            sth1(i) = tht1(i)
786
          enddo
787
      endif
3 michaesp 788
 
789
C       Determine the first and last loop indices
790
        if (numdat.eq.2) then
791
          filo = nint(tst/ts)+1
792
          lalo = nint((timeinc-ten)/ts) 
793
        elseif ( itm.eq.1 ) then
794
          filo = nint(tst/ts)+1
795
          lalo = nint(timeinc/ts)
796
        else if (itm.eq.numdat-1) then
797
          filo = 1
798
          lalo = nint((timeinc-ten)/ts)
799
        else
800
          filo = 1
801
          lalo = nint(timeinc/ts)
802
        endif
803
 
804
c       Split the interval <timeinc> into computational time steps <ts>
805
        do iloop=filo,lalo
806
 
807
C         Calculate relative time position in the interval timeinc (0=beginning, 1=end)
808
          reltpos0 = ((real(iloop)-1.)*ts)/timeinc
809
          reltpos1 = real(iloop)*ts/timeinc
810
 
811
c         Timestep for all trajectories
812
          do i=1,ntra
813
 
814
C           Check if trajectory has already left the data domain
815
            if (leftflag(i).ne.1) then	
816
 
19 michaesp 817
c             3D: Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
818
              if ( (imethod.eq.1   ).and.
819
     >             (isen   .eq.'no').and.
820
     >             (modlev .eq.'no') )
821
     >        then
822
                 call euler_3d(
3 michaesp 823
     >               xx1,yy1,pp1,leftflag(i),
824
     >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
825
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
826
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
827
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
828
 
19 michaesp 829
c             3D: Runge-Kutta timestep (x0,y0,p0 -> x1,y1,p1)
830
              else if ( (imethod.eq.2   ).and.
831
     >                  (isen   .eq.'no').and.
832
     >                  (modlev .eq.'no') )
833
     >        then
3 michaesp 834
                 call runge(
835
     >               xx1,yy1,pp1,leftflag(i),
836
     >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
837
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
838
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
839
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
840
 
19 michaesp 841
c             ISENTROPIC: Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
842
              else if ( (imethod.eq.1    ).and.
843
     >                  (isen   .eq.'yes').and.
844
     >                  (modlev .eq.'no' ) )
845
     >        then
846
                 call euler_isen(
847
     >               xx1,yy1,pp1,leftflag(i),
848
     >               xx0(i),yy0(i),pp0(i),theta(i),reltpos0,reltpos1,
849
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
850
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,
851
     >               sth0,sth1,tht0,tht1,
852
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
853
 
854
c             MODEL-LEVEL (2D): Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
855
              else if ( (imethod.eq.1    ).and.
856
     >                  (isen   .eq.'no' ).and.
857
     >                  (modlev .eq.'yes') )
858
     >        then
859
                 call euler_2d(
860
     >               xx1,yy1,pp1,leftflag(i),
861
     >               xx0(i),yy0(i),pp0(i),zindex(i),reltpos0,reltpos1,
862
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
863
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,
864
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
865
 
3 michaesp 866
              endif
867
 
868
c             Update trajectory position, or increase number of trajectories leaving domain
869
              if (leftflag(i).eq.1) then
870
                leftcount=leftcount+1
871
                if ( leftcount.lt.10 ) then
872
                   print*,'     -> Trajectory ',i,' leaves domain'
873
                elseif ( leftcount.eq.10 ) then
874
                   print*,'     -> N>=10 trajectories leave domain'
875
                endif
876
              else
877
                xx0(i)=xx1      
878
                yy0(i)=yy1
879
                pp0(i)=pp1
880
              endif
881
 
882
c          Trajectory has already left data domain (mark as <mdv>)
883
           else
884
              xx0(i)=mdv
885
              yy0(i)=mdv
886
              pp0(i)=mdv
887
 
888
           endif
889
 
890
          enddo
891
 
892
C         Save positions only every deltout minutes
893
          delta = aint(iloop*60*ts/deltout)-iloop*60*ts/deltout
894
          if (abs(delta).lt.eps) then
895
c          wstep = wstep + abs(ts)
896
c          if ( mod(wstep,deltout).eq.0 ) then
897
            time = time0+reltpos1*timeinc*fbflag
898
            itim = itim + 1
9 michaesp 899
            if ( itim.le.ntim ) then
900
              do i=1,ntra
901
                 call frac2hhmm(time,tload)
902
                 traout(i,itim,1) = tload
903
                 traout(i,itim,2) = xx0(i)
904
                 traout(i,itim,3) = yy0(i)
905
                 traout(i,itim,4) = pp0(i)
906
              enddo
907
            endif
3 michaesp 908
          endif
909
 
910
        enddo
911
 
912
      enddo
913
 
914
c     Write trajectory file
915
      vars(1)  ='time'
916
      vars(2)  ='lon'
917
      vars(3)  ='lat'
918
      vars(4)  ='p'
919
      call wopen_tra(cdfid,cdfname,ntra,ntim,4,reftime,vars,outmode)
920
      call write_tra(cdfid,traout,ntra,ntim,4,outmode)
921
      call close_tra(cdfid,outmode)   
922
 
923
c     Write some status information, and end of program message
924
      print*  
925
      print*,'---- STATUS INFORMATION --------------------------------'
926
      print*
927
      print*,'  #leaving domain    ', leftcount
928
      print*,'  #staying in domain ', ntra-leftcount
929
      print*
930
      print*,'              *** END OF PROGRAM CALTRA ***'
931
      print*,'========================================================='
932
 
933
      stop
934
 
935
c     ------------------------------------------------------------------
936
c     Exception handling
937
c     ------------------------------------------------------------------
938
 
939
 991  write(*,*) '*** ERROR: all start points outside the data domain'
940
      call exit(1)
941
 
942
 992  write(*,*) '*** ERROR: close arrays on files (prog. closear)'
943
      call exit(1)
944
 
945
 993  write(*,*) '*** ERROR: problems with array size'
946
      call exit(1)
947
 
948
      end 
949
 
950
 
951
c     *******************************************************************
952
c     * Time step : either Euler or Runge-Kutta                         *
953
c     *******************************************************************
954
 
955
C     Time-step from (x0,y0,p0) to (x1,y1,p1)
956
C
957
C     (x0,y0,p0) input	coordinates (long,lat,p) for starting point
958
C     (x1,y1,p1) output	coordinates (long,lat,p) for end point
959
C     deltat	 input	timestep in seconds
960
C     numit	 input	number of iterations
961
C     jump	 input  flag (=1 trajectories don't enter the ground)
962
C     left	 output	flag (=1 if trajectory leaves data domain)
963
 
964
c     -------------------------------------------------------------------
19 michaesp 965
c     Iterative Euler time step (KINEMATIC 3D TRAJECTORIES)
3 michaesp 966
c     -------------------------------------------------------------------
967
 
19 michaesp 968
      subroutine euler_3d(x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
969
     >                deltat,numit,jump,mdv,wfactor,fbflag,
970
     >		          spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
971
     >                xmin,ymin,dx,dy,per,hem,nx,ny,nz)
3 michaesp 972
 
973
      implicit none
974
 
975
c     Declaration of subroutine parameters
976
      integer      nx,ny,nz
977
      real         x1,y1,p1
978
      integer      left
979
      real	   x0,y0,p0
980
      real         reltpos0,reltpos1
981
      real   	   deltat
982
      integer      numit
983
      integer      jump
984
      real         wfactor
985
      integer      fbflag
986
      real     	   spt0(nx*ny)   ,spt1(nx*ny)
987
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
988
      real 	   vvt0(nx*ny*nz),vvt1(nx*ny*nz)
989
      real         wwt0(nx*ny*nz),wwt1(nx*ny*nz)
990
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
991
      real         xmin,ymin,dx,dy
992
      real         per
993
      integer      hem
994
      real         mdv
995
 
996
c     Numerical and physical constants
997
      real         deltay
998
      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
999
      real         pi                       
1000
      parameter    (pi=3.1415927)    ! Pi
1001
 
1002
c     Auxiliary variables
1003
      real         xmax,ymax
1004
      real	   xind,yind,pind
1005
      real	   u0,v0,w0,u1,v1,w1,u,v,w,sp
1006
      integer	   icount
1007
      character    ch
1008
 
1009
c     Externals    
1010
      real         int_index4
1011
      external     int_index4
1012
 
1013
c     Reset the flag for domain-leaving
1014
      left=0
1015
 
1016
c     Set the esat-north bounray of the domain
1017
      xmax = xmin+real(nx-1)*dx
1018
      ymax = ymin+real(ny-1)*dy
1019
 
1020
C     Interpolate wind fields to starting position (x0,y0,p0)
1021
      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
1022
     >                 p3d0,p3d1,spt0,spt1,3,
1023
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
1024
      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
1025
      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
1026
      w0 = int_index4(wwt0,wwt1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
1027
 
1028
c     Force the near-surface wind to zero
1029
      if (pind.lt.1.) w0=w0*pind
1030
 
1031
C     For first iteration take ending position equal to starting position
1032
      x1=x0
1033
      y1=y0
1034
      p1=p0
1035
 
1036
C     Iterative calculation of new position
1037
      do icount=1,numit
1038
 
1039
C        Calculate new winds for advection
1040
         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
1041
     >                    p3d0,p3d1,spt0,spt1,3,
1042
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
1043
         u1 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
1044
         v1 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
1045
         w1 = int_index4(wwt0,wwt1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
1046
 
1047
c        Force the near-surface wind to zero
1048
         if (pind.lt.1.) w1=w1*pind
1049
 
1050
c        Get the new velocity in between
1051
         u=(u0+u1)/2.
1052
         v=(v0+v1)/2.
1053
         w=(w0+w1)/2.
1054
 
1055
C        Calculate new positions
1056
         x1 = x0 + fbflag*u*deltat/(deltay*cos(y0*pi/180.))
1057
         y1 = y0 + fbflag*v*deltat/deltay
1058
         p1 = p0 + fbflag*wfactor*w*deltat/100.
1059
 
1060
c       Handle pole problems (crossing and near pole trajectory)
1061
        if ((hem.eq.1).and.(y1.gt.90.)) then
1062
          y1=180.-y1
1063
          x1=x1+per/2.
1064
        endif
1065
        if ((hem.eq.1).and.(y1.lt.-90.)) then
1066
          y1=-180.-y1
1067
          x1=x1+per/2.
1068
        endif
1069
        if (y1.gt.89.99) then
1070
           y1=89.99
1071
        endif
1072
 
1073
c       Handle crossings of the dateline
1074
        if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
1075
           x1=xmin+amod(x1-xmin,per)
1076
        endif
1077
        if ((hem.eq.1).and.(x1.lt.xmin)) then
1078
           x1=xmin+per+amod(x1-xmin,per)
1079
        endif
1080
 
1081
C       Interpolate surface pressure to actual position
1082
        call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
1083
     >                   p3d0,p3d1,spt0,spt1,3,
1084
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
1085
        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
1086
 
1087
c       Handle trajectories which cross the lower boundary (jump flag)
1088
        if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
1089
 
1090
C       Check if trajectory leaves data domain
1091
        if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
1092
     >       ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
1093
     >         (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
1094
     >  then
1095
          left=1
1096
          goto 100
1097
        endif
1098
 
1099
      enddo
1100
 
1101
c     Exit point for subroutine
1102
 100  continue
1103
 
1104
      return
1105
 
1106
      end
1107
 
1108
c     -------------------------------------------------------------------
19 michaesp 1109
c     Iterative Euler time step (ISENTROPIC)
1110
c     -------------------------------------------------------------------
1111
 
1112
      subroutine euler_isen(x1,y1,p1,left,x0,y0,p0,theta,
1113
     >                reltpos0,reltpos1,
1114
     >                deltat,numit,jump,mdv,wfactor,fbflag,
1115
     >                spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
1116
     >                sth0,sth1,tht0,tht1,
1117
     >                xmin,ymin,dx,dy,per,hem,nx,ny,nz)
1118
 
1119
      implicit none
1120
 
1121
c     Declaration of subroutine parameters
1122
      integer      nx,ny,nz
1123
      real         x1,y1,p1
1124
      integer      left
1125
      real         x0,y0,p0
1126
      real         reltpos0,reltpos1
1127
      real         deltat
1128
      integer      numit
1129
      integer      jump
1130
      real         wfactor
1131
      integer      fbflag
1132
      real         spt0(nx*ny)   ,spt1(nx*ny)
1133
      real         sth0(nx*ny)   ,sth1(nx*ny)
1134
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
1135
      real         vvt0(nx*ny*nz),vvt1(nx*ny*nz)
1136
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
1137
      real         tht0(nx*ny*nz),tht1(nx*ny*nz)
1138
      real         xmin,ymin,dx,dy
1139
      real         per
1140
      integer      hem
1141
      real         mdv
1142
      real         theta
1143
 
1144
c     Numerical and physical constants
1145
      real         deltay
1146
      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
1147
      real         pi
1148
      parameter    (pi=3.1415927)    ! Pi
1149
 
1150
c     Auxiliary variables
1151
      real         xmax,ymax
1152
      real         xind,yind,pind
1153
      real         u0,v0,w0,u1,v1,w1,u,v,w,sp
1154
      integer      icount
1155
      character    ch
1156
 
1157
c     Externals
1158
      real         int_index4
1159
      external     int_index4
1160
 
1161
c     Reset the flag for domain-leaving
1162
      left=0
1163
 
1164
c     Set the esat-north bounray of the domain
1165
      xmax = xmin+real(nx-1)*dx
1166
      ymax = ymin+real(ny-1)*dy
1167
 
1168
C     Interpolate wind fields to starting position (x0,y0,p0)
1169
      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
1170
     >                 p3d0,p3d1,spt0,spt1,3,
1171
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
1172
      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
1173
      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
1174
 
1175
C     For first iteration take ending position equal to starting position
1176
      x1=x0
1177
      y1=y0
1178
      p1=p0
1179
 
1180
C     Iterative calculation of new position
1181
      do icount=1,numit
1182
 
1183
C        Calculate new winds for advection
1184
         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
1185
     >                    p3d0,p3d1,spt0,spt1,3,
1186
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
1187
         u1 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
1188
         v1 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
1189
 
1190
c        Get the new velocity in between
1191
         u=(u0+u1)/2.
1192
         v=(v0+v1)/2.
1193
 
1194
C        Calculate new positions
1195
         x1 = x0 + fbflag*u*deltat/(deltay*cos(y0*pi/180.))
1196
         y1 = y0 + fbflag*v*deltat/deltay
1197
 
1198
c        Get the pressure on the isentropic surface at the new position
1199
         call get_index4 (xind,yind,pind,x1,y1,theta,reltpos1,
1200
     >                    tht0,tht1,sth0,sth1,1,
1201
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
1202
         p1 = int_index4(p3d0,p3d1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
1203
 
1204
c       Handle pole problems (crossing and near pole trajectory)
1205
        if ((hem.eq.1).and.(y1.gt.90.)) then
1206
          y1=180.-y1
1207
          x1=x1+per/2.
1208
        endif
1209
        if ((hem.eq.1).and.(y1.lt.-90.)) then
1210
          y1=-180.-y1
1211
          x1=x1+per/2.
1212
        endif
1213
        if (y1.gt.89.99) then
1214
           y1=89.99
1215
        endif
1216
 
1217
c       Handle crossings of the dateline
1218
        if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
1219
           x1=xmin+amod(x1-xmin,per)
1220
        endif
1221
        if ((hem.eq.1).and.(x1.lt.xmin)) then
1222
           x1=xmin+per+amod(x1-xmin,per)
1223
        endif
1224
 
1225
C       Interpolate surface pressure to actual position
1226
        call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
1227
     >                   p3d0,p3d1,spt0,spt1,3,
1228
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
1229
        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
1230
 
1231
c       Handle trajectories which cross the lower boundary (jump flag)
1232
        if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
1233
 
1234
C       Check if trajectory leaves data domain
1235
        if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
1236
     >       ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
1237
     >         (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
1238
     >  then
1239
          left=1
1240
          goto 100
1241
        endif
1242
 
1243
      enddo
1244
 
1245
c     Exit point for subroutine
1246
 100  continue
1247
 
1248
      return
1249
 
1250
      end
1251
 
1252
c     -------------------------------------------------------------------
1253
c     Iterative Euler time step (MODEL-LEVEL, 2D)
1254
c     -------------------------------------------------------------------
1255
 
1256
      subroutine euler_2d(x1,y1,p1,left,x0,y0,p0,zindex,
1257
     >                reltpos0,reltpos1,
1258
     >                deltat,numit,jump,mdv,wfactor,fbflag,
1259
     >                spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
1260
     >                xmin,ymin,dx,dy,per,hem,nx,ny,nz)
1261
 
1262
      implicit none
1263
 
1264
c     Declaration of subroutine parameters
1265
      integer      nx,ny,nz
1266
      real         x1,y1,p1
1267
      integer      left
1268
      real         x0,y0,p0
1269
      real         reltpos0,reltpos1
1270
      real         deltat
1271
      integer      numit
1272
      integer      jump
1273
      real         wfactor
1274
      integer      fbflag
1275
      real         spt0(nx*ny)   ,spt1(nx*ny)
1276
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
1277
      real         vvt0(nx*ny*nz),vvt1(nx*ny*nz)
1278
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
1279
      real         xmin,ymin,dx,dy
1280
      real         per
1281
      integer      hem
1282
      real         mdv
1283
      real         zindex
1284
 
1285
c     Numerical and physical constants
1286
      real         deltay
1287
      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
1288
      real         pi
1289
      parameter    (pi=3.1415927)    ! Pi
1290
      real         eps
1291
      parameter    (eps=0.001)
1292
 
1293
c     Auxiliary variables
1294
      real         xmax,ymax
1295
      real         xind,yind,pind
1296
      real         u0,v0,w0,u1,v1,w1,u,v,w,sp
1297
      integer      icount
1298
      character    ch
1299
 
1300
c     Externals
1301
      real         int_index4
1302
      external     int_index4
1303
 
1304
c     Reset the flag for domain-leaving
1305
      left=0
1306
 
1307
c     Set the esat-north bounray of the domain
1308
      xmax = xmin+real(nx-1)*dx
1309
      ymax = ymin+real(ny-1)*dy
1310
 
1311
C     Interpolate wind fields to starting position (x0,y0,p0)
1312
      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
1313
     >                 p3d0,p3d1,spt0,spt1,3,
1314
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
1315
      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,zindex,reltpos0,mdv)
1316
      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,zindex,reltpos0,mdv)
1317
 
1318
C     For first iteration take ending position equal to starting position
1319
      x1=x0
1320
      y1=y0
1321
      p1=p0
1322
 
1323
C     Iterative calculation of new position
1324
      do icount=1,numit
1325
 
1326
C        Calculate new winds for advection
1327
         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
1328
     >                    p3d0,p3d1,spt0,spt1,3,
1329
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
1330
         u1=int_index4(uut0,uut1,nx,ny,nz,xind,yind,zindex,reltpos1,mdv)
1331
         v1=int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,zindex,reltpos1,mdv)
1332
         if ( abs(u1-mdv).lt.eps ) then
1333
             left = 1
1334
             goto 100
1335
         endif
1336
         if ( abs(v1-mdv).lt.eps ) then
1337
             left = 1
1338
             goto 100
1339
         endif
1340
 
1341
c        Get the new velocity in between
1342
         u=(u0+u1)/2.
1343
         v=(v0+v1)/2.
1344
 
1345
C        Calculate new positions
1346
         x1 = x0 + fbflag*u*deltat/(deltay*cos(y0*pi/180.))
1347
         y1 = y0 + fbflag*v*deltat/deltay
1348
 
1349
c        Get the pressure on the model surface at the new position
1350
         xind = (x1 - xmin ) / dx + 1.
1351
         yind = (y1 - ymin ) / dy + 1.
1352
         p1 =
1353
     >     int_index4(p3d0,p3d1,nx,ny,nz,xind,yind,zindex,reltpos1,mdv)
1354
         if ( abs(p1-mdv).lt.eps ) then
1355
             left = 1
1356
             goto 100
1357
         endif
1358
 
1359
c       Handle pole problems (crossing and near pole trajectory)
1360
        if ((hem.eq.1).and.(y1.gt.90.)) then
1361
          y1=180.-y1
1362
          x1=x1+per/2.
1363
        endif
1364
        if ((hem.eq.1).and.(y1.lt.-90.)) then
1365
          y1=-180.-y1
1366
          x1=x1+per/2.
1367
        endif
1368
        if (y1.gt.89.99) then
1369
           y1=89.99
1370
        endif
1371
 
1372
c       Handle crossings of the dateline
1373
        if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
1374
           x1=xmin+amod(x1-xmin,per)
1375
        endif
1376
        if ((hem.eq.1).and.(x1.lt.xmin)) then
1377
           x1=xmin+per+amod(x1-xmin,per)
1378
        endif
1379
 
1380
C       Interpolate surface pressure to actual position
1381
        call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
1382
     >                   p3d0,p3d1,spt0,spt1,3,
1383
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
1384
        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
1385
 
1386
c       Handle trajectories which cross the lower boundary (jump flag)
1387
        if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
1388
 
1389
C       Check if trajectory leaves data domain
1390
        if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
1391
     >       ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
1392
     >         (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
1393
     >  then
1394
          left=1
1395
          goto 100
1396
        endif
1397
 
1398
      enddo
1399
 
1400
c     Exit point for subroutine
1401
 100  continue
1402
 
1403
      return
1404
 
1405
      end
1406
 
1407
c     -------------------------------------------------------------------
3 michaesp 1408
c     Runge-Kutta (4th order) time-step
1409
c     -------------------------------------------------------------------
1410
 
1411
      subroutine runge(x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
1412
     >                 deltat,numit,jump,mdv,wfactor,fbflag,
1413
     >		       spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
1414
     >                 xmin,ymin,dx,dy,per,hem,nx,ny,nz)
1415
 
1416
      implicit none
1417
 
1418
c     Declaration of subroutine parameters
1419
      integer      nx,ny,nz
1420
      real         x1,y1,p1
1421
      integer      left
1422
      real	   x0,y0,p0
1423
      real         reltpos0,reltpos1
1424
      real   	   deltat
1425
      integer      numit
1426
      integer      jump
1427
      real         wfactor
1428
      integer      fbflag
1429
      real     	   spt0(nx*ny)   ,spt1(nx*ny)
1430
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
1431
      real 	   vvt0(nx*ny*nz),vvt1(nx*ny*nz)
1432
      real         wwt0(nx*ny*nz),wwt1(nx*ny*nz)
1433
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
1434
      real         xmin,ymin,dx,dy
1435
      real         per
1436
      integer      hem
1437
      real         mdv
1438
 
1439
c     Numerical and physical constants
1440
      real         deltay
1441
      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
1442
      real         pi                       
1443
      parameter    (pi=3.1415927)    ! Pi
1444
 
1445
c     Auxiliary variables
1446
      real         xmax,ymax
1447
      real	   xind,yind,pind
1448
      real	   u0,v0,w0,u1,v1,w1,u,v,w,sp
1449
      integer	   icount,n
1450
      real         xs,ys,ps,xk(4),yk(4),pk(4)
1451
      real         reltpos
1452
 
1453
c     Externals    
1454
      real         int_index4
1455
      external     int_index4
1456
 
1457
c     Reset the flag for domain-leaving
1458
      left=0
1459
 
1460
c     Set the esat-north bounray of the domain
1461
      xmax = xmin+real(nx-1)*dx
1462
      ymax = ymin+real(ny-1)*dy
1463
 
1464
c     Apply the Runge Kutta scheme
1465
      do n=1,4
1466
 
1467
c       Get intermediate position and relative time
1468
        if (n.eq.1) then
1469
          xs=0.
1470
          ys=0.
1471
          ps=0.
1472
          reltpos=reltpos0
1473
        else if (n.eq.4) then
1474
          xs=xk(3)
1475
          ys=yk(3)
1476
          ps=pk(3)
1477
          reltpos=reltpos1
1478
        else
1479
          xs=xk(n-1)/2.
1480
          ys=yk(n-1)/2.
1481
          ps=pk(n-1)/2.
1482
          reltpos=(reltpos0+reltpos1)/2.
1483
        endif
1484
 
1485
C       Calculate new winds for advection
1486
        call get_index4 (xind,yind,pind,x0+xs,y0+ys,p0+ps,reltpos,
1487
     >                   p3d0,p3d1,spt0,spt1,3,
1488
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
1489
        u = int_index4 (uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos,mdv)
1490
        v = int_index4 (vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos,mdv)
1491
        w = int_index4 (wwt0,wwt1,nx,ny,nz,xind,yind,pind,reltpos,mdv)
1492
 
1493
c       Force the near-surface wind to zero
1494
        if (pind.lt.1.) w1=w1*pind
1495
 
1496
c       Update position and keep them
1497
        xk(n)=fbflag*u*deltat/(deltay*cos(y0*pi/180.))
1498
        yk(n)=fbflag*v*deltat/deltay
1499
        pk(n)=fbflag*w*deltat*wfactor/100.
1500
 
1501
      enddo
1502
 
1503
C     Calculate new positions
1504
      x1=x0+(1./6.)*(xk(1)+2.*xk(2)+2.*xk(3)+xk(4))
1505
      y1=y0+(1./6.)*(yk(1)+2.*yk(2)+2.*yk(3)+yk(4))
1506
      p1=p0+(1./6.)*(pk(1)+2.*pk(2)+2.*pk(3)+pk(4))
1507
 
1508
c     Handle pole problems (crossing and near pole trajectory)
1509
      if ((hem.eq.1).and.(y1.gt.90.)) then
1510
         y1=180.-y1
1511
         x1=x1+per/2.
1512
      endif
1513
      if ((hem.eq.1).and.(y1.lt.-90.)) then
1514
         y1=-180.-y1
1515
         x1=x1+per/2.
1516
      endif
1517
      if (y1.gt.89.99) then
1518
         y1=89.99
1519
      endif
1520
 
1521
c     Handle crossings of the dateline
1522
      if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
1523
         x1=xmin+amod(x1-xmin,per)
1524
      endif
1525
      if ((hem.eq.1).and.(x1.lt.xmin)) then
1526
         x1=xmin+per+amod(x1-xmin,per)
1527
      endif
1528
 
1529
C     Interpolate surface pressure to actual position
1530
      call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
1531
     >                 p3d0,p3d1,spt0,spt1,3,
1532
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
1533
      sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1,reltpos,mdv)
1534
 
1535
c     Handle trajectories which cross the lower boundary (jump flag)
1536
      if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
1537
 
1538
C     Check if trajectory leaves data domain
1539
      if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
1540
     >     ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
1541
     >       (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
1542
     >then
1543
         left=1
1544
         goto 100
1545
      endif
1546
 
1547
c     Exit point fdor subroutine
1548
 100  continue
1549
 
1550
      return
1551
      end