Subversion Repositories lagranto.icon

Rev

Rev 3 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
 
2
      SUBROUTINE special_stt (flag,cmd,tra,ntim,ncol,
3
     >                    vars,times,param,nparam,outfile)
4
 
5
c     ***************************************************************************
6
c     *                                                                         *
7
c     * OUTPUT:  flag           -> 1 if trajectory is selected, 0 if not        *
8
c     *                                                                         *
9
c     * INPUT:   cmd            <- command string                               *
10
c     *          tra(ntim,ncol) <- single trajectory: indices time,column       *
11
c     *          ntim           <- number of times                              *
12
c     *          ncol           <- number of columns (including time,lon,lat,p) *
13
c     *          vars(ncol)     <- names of columns                             *
14
c     *          times(ntim)    <- List of times
15
c     *          param(nparam)  <- parameter values                             *
16
c     *          nparam         <- number of parameters                         *
17
c     *                                                                         *
18
c     ***************************************************************************
19
 
20
      implicit none
21
 
22
c     ---------------------------------------------------------------------------
23
c     Declaration of subroutine parameters
24
c     ---------------------------------------------------------------------------
25
 
26
      integer       flag           ! Boolean flag whether trajectory is selected
27
      character*80  cmd            ! Command string
28
      integer       ntim,ncol      ! Dimension of single trajectory
29
      real          tra(ntim,ncol) ! Single trajectory
30
      character*80  vars(ncol)     ! Name of columns
31
      real          times(ntim)    ! List of times
32
      integer       nparam         ! # parameters
33
      real          param(1000)    ! List of parameters
34
      character*80  outfile        ! NAme of output file
35
 
36
      integer       i,n		   ! Loop indices
37
 
38
c     ---------------------------------------------------------------------------
39
c     Local variables for SPECIAL:WCB
40
c     ---------------------------------------------------------------------------
41
 
42
      integer       ip,i0,i1
43
 
44
c     ---------------------------------------------------------------------------
45
c     Local variables for SPECIAL:STT
46
c     ---------------------------------------------------------------------------
47
 
48
      integer       debugflag
49
      integer       tb,ta					!tb and ta are the minimum residence times before and after the exchange
50
      integer	    stat, ex, outputflag			!Variables for I/O files
51
      integer	    labelflag					!Flag to handle label field
52
      integer	    pvflag, thflag, zflag, pflag		!Flags to signal presence of various fields
53
      integer	    pvcol,thcol,labelcol,zcol,pcol		!Positions of PV, TH, LABEL, z and p columns
54
      real	    mean1, mean2, sd1, sd2			!Mean and standard deviations of PV and TH
55
      real	    deltat					!Trajectory timestep
56
      real	    extime, exlon, exlat, exth, exheight, expres!Position at exchange time	
57
      character*80  namefile            			!Name of output file
58
      real	    a,b						!Normalized tropopause distance for PV and TH
59
      real	    thtrop, pvtrop				!Tropopause level in K and pvu
60
      integer	    count1, count2
61
      integer	    sttcount, p				!Ancillary variables in STT do while loop
62
      real,allocatable, dimension (:)   :: tropdist	! Vector for normalized distance from tropopause
63
      real,allocatable, dimension (:)   :: exflag	! Vector for normalized distance from tropopause
64
      integer, dimension (ntim)		:: pos		! Vector to store exchange indices
65
 
66
 
67
 
68
c     --------------------------------------------------------------------------  %)
69
c     SPECIAL:STT:pvtrop,thtrop,tmin_b,tmin_f,outputflag 		                  	  %)
70
 
71
c     --------------------------------------------------------------------------- %)
72
 
73
      if ( cmd.eq.'STT' ) then
74
 
75
c        Reset the flag for selection
76
         flag = 0
77
c	 Reset the flag for debugging
78
	 debugflag = 0
79
c	 Set other flags
80
         labelflag = 0
81
	 pvflag = 0
82
	 thflag = 0
83
C	 Set data at the exchange
84
         extime = 0
85
	 exlon = 0 
86
	 exlat = 0
87
	 exheight = 0
88
	 expres = 0
89
 
90
c	 PV and TH threshold read from command line
91
	 pvtrop = param(1)
92
	 thtrop = param(2)
93
 
94
 
95
c	Determine number of time steps in the trajectory 
96
c	Handles backward OR forward trajectories
97
 
98
	if ( (ntim.eq.0).or.(ntim.eq.1) ) then
99
          print*,' ERROR: zero length trajectory in SPECIAL:STT... Stop'	!If there are problems, stop computation.
100
            stop
101
         endif
102
 
103
	deltat = (times(ntim)-times(1))/(ntim-1)				!By default array indexing starts in Fortran at 1
104
 
105
 
106
C 	Now some checks for essential fields : label, potential vorticity, potential temperature
107
 
108
c	Check for the presence of the LABEL field
109
	labelcol = 1
110
	do i = 1,ncol
111
	    if((vars(i) .eq. 'LABEL').or.(vars(i) .eq. 'label')) then
112
		labelflag = 1
113
		labelcol = i
114
	    if(debugflag .eq. 1) then
115
		print*, 'Found label'
116
	    endif
117
	    endif
118
	enddo
119
 
120
 
121
c	Check for the presence of the variables PV and TH
122
 
123
	pvcol = 1	
124
	thcol = 1
125
	zcol = 1
126
	pcol = 1			!In case th or pv are not available, tra(timestep,thcol) is always defined
127
 
128
 
129
	do i = 1,ncol
130
	    if((vars(i) .eq. 'PV') .or. (vars(i) .eq. 'POT_VORTI')) then
131
		pvflag = 1
132
		pvcol = i
133
	    endif
134
	enddo
135
 
136
 
137
	do i = 1,ncol
138
	    if((vars(i) .eq. 'TH') .or. (vars(i) .eq. 'POT_T')) then
139
		thflag = 1
140
		thcol = i
141
	    endif
142
	enddo
143
 
144
	do i = 1,ncol
145
	    if((vars(i) .eq. 'z') .or. (vars(i) .eq. 'Z')) then
146
		zcol = i
147
		zflag = 1
148
	    endif
149
	enddo
150
 
151
	do i = 1,ncol
152
	    if((vars(i) .eq. 'p') .or. (vars(i) .eq. 'P')) then
153
		pcol = i
154
		pflag = 1
155
	    endif
156
	enddo
157
 
158
c	No PV and THETA available: not possible determining tropopause. Program stop. 
159
 
160
	if ( (thflag .eq. 0).and.(pvflag .eq. 0) ) then
161
            print*,' ERROR: no PV and no TH in SPECIAL:STT..'
162
	    print*,' ERROR: not possible determining tropopause'	
163
            stop
164
         endif
165
 
166
c	Check for missing data in the trajectories
167
 
168
  	do n=1,ntim
169
     		if ((tra(n,pvcol).lt.-500.).and.(pvflag.eq.1.)) then
170
	          if(debugflag .eq. 1) then
171
            		print*,'WARNING: PV missing data',tra(n,pvcol)
172
			print*,'Trajectory skipped'
173
		  endif
174
			goto 600
175
    		 endif 
176
		if ((tra(n,thcol) .lt.-500.).and.(thflag .eq. 1.)) then
177
		  if(debugflag .eq. 1) then
178
            		print*,'WARNING: TH missing data',tra(n,thcol)
179
			print*,'Trajectory skipped'
180
		  endif
181
			goto 600
182
    		 endif 
183
  	enddo
184
 
185
c	Computing mean of PV and TH along a trajectory
186
 
187
	mean1=0.
188
  	count1=0
189
  	mean2=0.
190
  	count2=0
191
  	do n=1,ntim						
192
     	    if ((tra(n,pvcol).gt.-100.).and.(pvflag .eq. 1)) then
193
        	count1=count1+1
194
        	mean1=mean1+tra(n,pvcol)
195
	    else
196
	    if(debugflag .eq. 1) then
197
		print*,'WARNING: PV is missing or non-physical, mean1 = 0'
198
	    endif
199
		mean1 = 0
200
	    endif
201
     	    if ((tra(n,thcol) .gt. 1E-6) .and. (thflag .eq. 1)) then
202
        	count2=count2+1
203
        	mean2=mean2+tra(n,thcol)
204
	    else
205
	    if(debugflag .eq. 1) then
206
		print*,'WARNING: TH is missing or non-physical, mean1 = 0'
207
	    endif
208
		mean2 = 0
209
	    endif
210
  	enddo
211
 
212
c	Compute average, otherwise skip problematic trajectories
213
 
214
  	if (count1.ne.0) then
215
     		mean1=mean1/count1
216
  	endif
217
 
218
  	if (count2.ne.0) then
219
     		mean2=mean2/count2
220
  	endif
221
 
222
        if ( (count1.eq.0) .and. (count2.eq.0) )then
223
	    if(debugflag .eq. 1) then
224
		print*, 'Trajectory with only missing data skipped' 
225
	    endif       	
226
		goto 600
227
        endif
228
 
229
c	Computing standard deviation of PV and TH along a trajectory
230
 
231
	sd1=0.
232
  	count1=0.
233
  	sd2=0.
234
  	count2=0.
235
  	do n=1,ntim						
236
     	    if (tra(n,pvcol).gt.-100.) then
237
        	count1=count1+1.
238
        	sd1=sd1+(tra(n,pvcol)-mean1)**2
239
	    endif
240
     	    if (tra(n,thcol).gt. 0.) then
241
        	count2=count2+1.
242
        	sd2=sd2+(tra(n,thcol)-mean2)**2
243
	    endif
244
  	enddo
245
	sd1=sqrt(sd1/(count1-1.))
246
  	sd2=sqrt(sd2/(count2-1.))
247
 
248
c     Allocate distance array
249
 
250
	allocate(tropdist(1:ntim),stat=stat)
251
	if (stat.ne.0) then 
252
	    print*,'ERROR: allocating array tropdist 
253
     &not possible in SPECIAL:STT.. '
254
	    stop
255
	endif
256
 
257
c	If LABEL is not available, go to section without usage of label
258
 
259
	if(labelflag .ne. 1) goto 300
260
 
261
 
262
C-------------------------------------------------------------------------------------------------------------------------
263
c	BEGINNING OF THE SECTION WITH LABEL
264
C-------------------------------------------------------------------------------------------------------------------------
265
 
266
 
267
c	Check for label field
268
	do i = 1,ntim
269
	    if((tra(i,labelcol) .ne.1) .and. (tra(i,labelcol) .ne.2)) then
270
	      if(debugflag .eq. 1) then
271
	    	print*,'WARNING: unexpected label field in SPECIAL:STT.. '
272
	      endif
273
		print*,'WARNING: unused label field in SPECIAL:STT'
274
	    	goto 300
275
	    endif
276
	end do
277
 
278
 
279
c     Distance from the '2 PVU, 380 K' tropopause (-:trop/+:strat) : pvtrop, thtrop if PV and TH and LABEL are available
280
 
281
	if ( (pvflag .ne. 0) .and. (thflag .ne. 0) .and. 
282
     &(labelflag .eq. 1)) then 
283
	    if(debugflag .eq. 1) then
284
	      print*,'Tropopause detected with PV, TH and LABEL'
285
            endif
286
	  do n=1,ntim
287
     	    a=abs(pvtrop-abs(tra(n,pvcol)))/sd1 	! a = tropdist(PV)
288
     	    b=abs(thtrop-tra(n,thcol))/sd2 		! b = tropdist(TH)
289
     	    if (tra(n,labelcol) .eq. 1) then ! Troposphere (LABEL)  
290
        	tropdist(n)=-min(a,b)
291
     	    else
292
         	tropdist(n)=+min(a,b)
293
     	    endif
294
  	  enddo
295
	endif
296
 
297
c	If PV is not available, use only TH
298
 
299
	if  (pvflag .eq. 0) then 
300
	  if(debugflag .eq. 1) then
301
	    print*,'Tropopause detected with TH only'
302
	  endif
303
	  do n=1,ntim
304
     	    b=abs(thtrop-tra(n,thcol))/sd2 		! b = tropdist(TH)
305
     	    if(tra(n,labelcol) .eq. 1) then 		! Troposphere (TH)
306
        	tropdist(n)= -b
307
	    else
308
	        tropdist(n)= b				! Stratosphere (TH)
309
     	    endif
310
  	  enddo
311
	endif
312
 
313
c	If TH is not available, use only PV
314
	if  (thflag .eq. 0) then
315
	  if(debugflag .eq. 1) then
316
	    print*,'Tropopause detected with PV only'	
317
	  endif
318
	  do n=1,ntim
319
     	    a=abs(pvtrop-tra(n,pvcol))/sd2 		! a = tropdist(PV)
320
     	    if( (tra(n,pvcol).lt.pvtrop) .and. 
321
     &(tra(n,labelcol) .eq. 1)) then 		! Troposphere (PV)
322
        	tropdist(n)= -a
323
     	    else
324
	        tropdist(n)= a				! Stratosphere (PV)
325
     	    endif
326
  	  enddo
327
	endif
328
 
329
C	Detecting STT with LABEL field
330
c	If PV is lower than pvtrop despite the label is 2, then in next step this potential STT is not considered.
331
 
332
	sttcount = 0
333
 
334
	do i=2,ntim
335
	    if((tropdist(i)*tropdist(i-1) .lt. 0) .and. 
336
     &(tropdist(i-1) .gt. 0)) then
337
C		Adding label check for STT
5 michaesp 338
		if((tra(i,labelcol).ne.1).and.(tra(i-1,labelcol).ne.2))then
3 michaesp 339
		  if(debugflag .eq. 1) then
340
	            print*, 'STT label check failed:' 
341
		    print*, 'moving to next STT candidate'
342
	          endif
343
		    cycle
344
		else
345
		    sttcount = sttcount+1
346
		    pos(sttcount) = i-1
347
		endif
348
	    endif
349
	enddo
350
 
351
c	If no STT is detected, move to next trajectory
352
 
353
	if(sttcount .lt. 1 ) then
354
	  if(debugflag .eq. 1) then
355
	    print*, 'No STT detected: moving to next trajectory'
356
	  endif
357
	    goto 600
358
	endif
359
 
360
c	Re-connection with normal steps: label has been used
361
 
362
	goto 350
363
 
364
C-------------------------------------------------------------------------------------------------------------------------
365
C	END OF LABEL SECTION
366
c-------------------------------------------------------------------------------------------------------------------------
367
 
368
 
369
 
370
c     Distance from the '2 PVU, 380 K' tropopause (-:trop/+:strat) : pvtrop, thtrop if PV and TH are available
371
300	continue							!Here the code continues if no LABEL has been found
372
 
373
	if ( (pvflag .ne. 0) .and. (thflag .ne. 0)) then
374
	    if(debugflag .eq. 1) then
375
	      print*,'Tropopause detected with PV and TH'
376
            endif
377
	  do n=1,ntim
378
     	    a=abs(pvtrop-abs(tra(n,pvcol)))/sd1 	! a = tropdist(PV)
379
     	    b=abs(thtrop-tra(n,thcol))/sd2 		! b = tropdist(TH)
380
     	    if ((tra(n,thcol).lt.thtrop).and.
381
     &(abs(tra(n,pvcol)).lt.pvtrop)) then ! Troposphere (TH and PV)	 ---------WHY ABS??---------
382
        	tropdist(n)=-min(a,b)
383
     	    else
384
       	    	if ( (abs(tra(n,pvcol)).ge.pvtrop).and.
385
     &(tra(n,thcol).ge.thtrop) ) then ! Stratosphere (TH and PV)
386
           	    tropdist(n)=+min(a,b)
387
            	endif
388
 
389
            	if ( (abs(tra(n,pvcol)).ge.pvtrop).and.
390
     &(tra(n,thcol).lt.thtrop) ) then ! Stratosphere (PV)
391
           	    tropdist(n)=a
392
            	endif
393
 
394
            	if ( (abs(tra(n,pvcol)).lt.pvtrop).and.
395
     &(tra(n,thcol).ge.thtrop) ) then ! Stratosphere (TH)
396
           	    tropdist(n)=b
397
            	endif
398
     	    endif
399
  	  enddo
400
	endif
401
 
402
c	If PV is not available, use only TH
403
 
404
	if  (pvflag .eq. 0) then
405
	    if(debugflag .eq. 1) then
406
	      print*,'Tropopause detected with TH only'
407
            endif
408
	  do n=1,ntim
409
     	    b=abs(thtrop-tra(n,thcol))/sd2 		! b = tropdist(TH)
410
     	    if (tra(n,thcol).lt.thtrop) then 		! Troposphere (TH)
411
        	tropdist(n)= -b
412
     	    else
413
	        tropdist(n)= b				! Stratosphere (TH)
414
     	    endif
415
  	  enddo
416
	endif
417
 
418
c	If TH is not available, use only PV
419
	    if(debugflag .eq. 1) then
420
	      print*,'Tropopause detected with PV only'
421
            endif	
422
	if  (thflag .eq. 0) then
423
	  do n=1,ntim
424
     	    a=abs(pvtrop-tra(n,pvcol))/sd2 		! a = tropdist(PV)
425
     	    if (tra(n,pvcol).lt.pvtrop) then 		! Troposphere (PV)
426
        	tropdist(n)= -a
427
     	    else
428
	        tropdist(n)= a				! Stratosphere (PV)
429
     	    endif
430
  	  enddo
431
	endif
432
 
433
 
434
c	Detect the various STT occurrences and save positions in traj file in array pos
435
 
436
	sttcount = 0
437
 
438
	do i=2,ntim
439
	    if((tropdist(i)*tropdist(i-1) .lt. 0) .and.
440
     & (tropdist(i-1) .gt. 0)) then
441
		sttcount = sttcount+1
442
		pos(sttcount) = i-1
443
	    endif
444
	enddo
445
 
446
c	If no STT are detected, go to next trajectory
447
 
448
	if(sttcount .lt. 1 ) then
449
	  if(debugflag .eq. 1) then
450
	    print*, 'No STT detected: moving to next trajectory'
451
	  endif
452
	    goto 600
453
	endif
454
 
455
 
456
 
457
350	continue								!The same if LABEL is present or not.
458
 
459
c	In timesteps, minimum residence time before and after the exchange
460
 
461
	tb = nint( param(3)/deltat )
462
	ta = nint( param(4)/deltat )
463
 
464
 
465
c	Now a big while loop analyzing every candidate STT
466
 
467
	do while(sttcount>0)
468
 
469
	    p = pos(sttcount)	
470
 
471
c	    Check whether the trajectory is still in the troposphere after ta timesteps
472
 
473
	   if(ta .ne. 0) then
474
	    	do n=1,ta
475
	            if((p+n) .le. ntim) then
476
		    	if (tropdist(p)*tropdist(p+n).gt.0.) goto 500		!Try with new STT candidate
477
		    else
478
		    	goto 500						!Out of range. Try with new STT candidate							
479
	    	    endif
480
  	    	enddo
481
	   endif
482
 
483
c	    Check whether the trajectory had been in the stratosphere for at least tb timesteps
484
 
485
	    if(tb .ne. 0) then
486
	        do n=1,tb
487
	            if((p-n) .gt. 0) then
488
		    	if (tropdist(p)*tropdist(p-n).lt.0.) goto 500		!Try with new STT candidate
489
		    else
490
		    	goto 500						!Out of range. Try with new STT candidate
491
	    	    endif
492
  	    	enddo
493
	    endif
494
 
495
c	   If the loops have gone fine, the trajectory is considered, flag = 1, values at the exchange are written and the loop stops.
496
 
497
	    flag = 1
498
	    extime = tra(p,1)
499
	    exlon = tra(p,2)
500
	    exlat = tra(p,3)
501
	    if(zflag .eq. 1) then
502
	    	exheight = tra(p,zcol)
503
	    endif
504
	    if(pflag .eq. 1) then
505
	    	expres = tra(p,pcol)
506
	    endif
507
	    if(thflag .eq. 1) then
508
	    	exth = tra(p,thcol)
509
	    endif
510
	    goto 600
511
 
512
500	    continue
513
	    sttcount = sttcount-1 
514
 
515
	enddo		! End of do while loop
516
 
517
	    goto 600
518
 
519
 
520
 
521
600   continue
522
 
523
      outputflag = param(5)
524
 
525
c 	Read the name of the file
526
 
527
	    namefile = trim(outfile)//'.stt'
528
 
529
 	if( (flag .eq. 1) .and. (outputflag .eq. 1) ) then
530
 
531
      if ( nint( param(nparam+1) ).eq.0) then
532
		open(15,file=namefile,status='replace',	
533
     &position='append', iostat=stat)
534
 
535
c        ---------------------------------------------------------------
536
       write(15,*) 'time      lon     lat     height  p       th'
537
       write(15,*) '--------------------------------------------'
538
       write(15,*) 
539
            	write(15,'(f7.2,1x,2f8.2,f8.0,9f10.3,5x,i1,4x)') 
540
     &extime,exlon,exlat,exheight,expres,exth
541
c        ---------------------------------------------------------------
542
 
543
        	close(15)
544
            param(nparam+1) = 1.
545
 
546
      else
547
 
548
		open(15,file=namefile,status='old',	
549
     &position='append', iostat=stat)
550
 
551
c        ---------------------------------------------------------------
552
            	write(15,'(f7.2,1x,2f8.2,f8.0,9f10.3,5x,i1,4x)') 
553
     &extime,exlon,exlat,exheight,expres,exth
554
c        ---------------------------------------------------------------
555
 
556
        	close(15)
557
      endif
558
 
559
 
560
 
561
 
562
 	endif 
563
 
564
 
565
c	Write output if outputflag = 1
566
 
567
c	outputflag = param(5)
568
 
569
c	Check for correct inserted parameter
570
 
571
c	if(debugflag .eq. 1) then
572
c	if( (outputflag .ne. 0) .and. (outputflag .ne. 1) ) then
573
c	  print*, 'WARNING: last parameter must be either 0 or 1'
574
c	  print*, 'NO data file produced'
575
c	endif
576
c	endif
577
 
578
c        if( (flag .eq. 1) .and. (outputflag .eq. 1) ) then
579
 
580
c	Read the name of the file
581
c	    namefile = "data_exch_stt"							!Maybe write in the command 'name'
582
c	Open file (already opened before)
583
c	
584
c            inquire(file=namefile, exist=ex)
585
c	    if (ex) then
586
c		open(15,file=namefile, status = 'old', 
587
c     &position='append',iostat=stat)
588
c	    if(stat .eq. -1) then							!If the file is empty behave as it was non-existent
589
c	      goto 125
590
c	    endif
591
c        Write output for trajectory making STT
592
c        ---------------------------------------------------------------
593
c            	write(15,'(f7.2,1x,2f8.2,f6.0,9f10.3,5x,i1,4x)') 
594
c     &extime,exlon,exlat,exheight,expres,exth
595
c        ---------------------------------------------------------------
596
c            	write(15,*) ! empty line
597
c
598
c
599
c	    else
600
c
601
c	First time opening
602
c
603
c		open(15,file=namefile,status='replace',
604
c     &position='append', iostat=stat)
605
c
606
c        Write title
607
c        ---------------------------------------------------------------
608
c 125           	write(15,"(A)") "time	lon	lat	z	p	th"
609
c          write(15,"(A)") "-------------------------------------------"
610
 
611
c        ---------------------------------------------------------------
612
c            	write(15,*) ! empty line
613
 
614
c        Write output for trajectory making STT
615
c        ---------------------------------------------------------------
616
c            	write(15,'(f7.2,1x,2f8.2,f6.0,9f10.3,5x,i1,4x)')
617
c     &extime,exlon,exlat,exheight,expres,exth
618
c        ---------------------------------------------------------------
619
c            	write(15,*) ! empty line
620
c	    endif
621
 
622
 
623
c	    if (stat .ne. 0) then
624
c		print *,"ERROR: Problem with output file",namefile
625
c	        stop
626
c	    endif
627
c	close(15)
628
c	endif
629
 
630
c	"Closing" the STT criterion
631
	endif 
632
 
633
       end