-
Notifications
You must be signed in to change notification settings - Fork 1
/
make_tracks.ftn
2551 lines (2335 loc) · 98.5 KB
/
make_tracks.ftn
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
program make_tracks
implicit none
#ifdef NETCDF
include 'netcdf.inc'
#endif
!
!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
! ----------------
! Make_Tracks v2.3
! ----------------
! Programme pour determiner de facon objective les trajectoires
! des depressions et anticyclones.
!
!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
! > Code porte sous SVN en Septembre 2015.
!
! **** VERSIONS PRE SVN ... ***
!
! > Version 2.2 - J.-F. Caron, Aout 2014
! -----------------------------------------
! - Mode WarmStart permettant de specifier les positions des
! cyclones actifs au temps initial
! - Option "-onlywarm" pour suivre seulement les cyclones
! specifies dans le fichier ./warmstart.txt
! - Prise en compte du masque Terre-Mer pour les cyclones tropicaux
!
! > Version 2.1 - J.-F. Caron, Septembre 2011
! -----------------------------------------
! - Ajout du suivi des cyclones tropicaux et des diagnostics pour
! determiner la phase des cyclones. Idees tirees des articles de
! M. R. Sinclair (2004, MWR) et de R. E. Hart (2003, MWR),
! ainsi que du programme track_hurricane.F (par K. Winger,
! L.-P. Caron et A.-M. Leduc).
! -----------------------------------------
!
! > Version 2.0 - J.-F. Caron, Juillet 2011
! -----------------------------------------
! - Refonte majeure du code : version non-implicite,
! allocation dynamique des vecteurs lies a la grille, utilisation
! des fonctions EZSCINT, support des grilles Lat-Lon
! incluant les grilles globales, redistribution du code entre 2
! librairies ( libtracks et libtracks_tools ).
! -----------------------------------------
!
! > Version 1.4 - J.-F. Caron, Automne 2010
! -----------------------------------------
! - Reformulation du calcul des probabilites (variable PROB).
! - Extension pour attenuer les effets du filtre sur le tourbillon
! en bordure de la grille (pour les LAM).
! - Imposition de valeurs seuils sur le changement de direction
! de deplacement.
! - Ajustements pour suivi en pression.
! -----------------------------------------
!
! > Version 1.3 - J.-F. Caron, Printemps-Ete 2007
! -----------------------------------------
! - Filtrage deplace APRES le calcul du tourbillon.
! - Re-organisation partielle du code.
! (e.g. regroupement des S-R dans la librairie LIBTRACKS)
! -----------------------------------------
!
! > Version 1.2 - Milka Radojevic, SCA/UQAM, Octobre 2004
! -----------------------------------------
! - Correction du calcul de la circulation.
! [ Activite des cyclones extra-tropicaux simules par le modele
! couple canadian de circulation globale (MCCG3)" - Memoire de
! maitrise, UQAM, 2006 ]
! -----------------------------------------
!
! > Version 1.1 - Corina Rosu Costea, SCA/UQAM, Nov 2003-Avr 2004
! -----------------------------------------
! - Adaptation au format RPN.
! - On prend la moitie du vent a 500 hPa au lieu du vent
! climatologique pour estimer la vitesse de deplacement.
! [ Les caracteristiques des cyclones et l'apport d'eau dans les
! bassins versants du Quebec" - Memoire de maitrise, UQAM, 2005 ]
! -----------------------------------------
!
! > Version 1.0 - Aout 2003
! -----------------------------------------
! - Reception du code de Sinclair
! [ Mark R. Sinclair, 1997: "Objective Identification of Cyclones
! and Their Circulation Intensity, and Climatology", Weather and
! Forecasting, 12, 3, 595-612. ]
! -----------------------------------------
!
!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
#include "tracks_cte.h"
!
integer maxtim, maxtrk, maxctr
parameter(maxtim=40000) ! Nombre de pas de temps max.
parameter(maxtrk=500) ! Nombre de trajectoires actives max.
parameter(maxctr=150) ! Nombre de centres max. par trajectoire
! Variables pour les OPTIONS
integer naj ! Nombre de pas de temps par jour (24 hrs)
integer npt ! Temps de vie minimum (en pas de temps)
! d'une trajectoire
integer level ! Niveau vertical suivi
integer fvort ! Cressman smoothing radius (km) for vort
integer fpres ! Cressman smoothing radius (km) for press
integer analyse ! >= 0 ( < 0) analyses (previsions) en entree
integer nofst ! /=0 ( =0 pas d') ecriture en fichier std ou NetCDF
integer isign ! +1 (-1) on suit les cyclones (anticyclones)
integer vortcnt ! 0 ( /=0) suivi en tourbillon (pression)
integer anticycl ! 0 ( /=0) on suit les anticyclones (cyclones)
integer phase ! >=0 on diagnotic la phase des cyclones
integer tropic ! >=0 on suit uniquement les cyclones tropicaux
integer latmin ! Frontiere Sud
integer latmax ! Frontiere Nord
integer lonmin ! Frontiere Ouest
integer lonmax ! Frontiere Est
integer relocp ! 0 ( /=0) on (ne) repositionne (pas) les
! centres de pression filtree
integer onlywarm ! 0 = On suit seulement les cyclones specifies
! dans le fichier ./warmstart.txt
integer silent ! Regle le niveau d'ecriture a l'ecran
integer hcrit ! Valeur seuil pour les epaisseurs
integer bcrit ! Valeur seuil pour le terme de baroclinicite
integer tcrad ! Rayon d'un cyclone tropical
integer tdwind ! Valeur seuil de vent pour considerer un
! cyclone comme depression tropicale
integer slevel ! Niveau pour le vent de surface (10m)
real dist_crit ! Rayon du cercle de recherche (en deg)
real vc ! Valeur seuil pour le tourbillon
logical vort_c ! True (False) : on suit les extremes
! de tourbillon (pression)
logical do_circ ! True : on calcule la circulation
character*12 output_type ! format des sorties ascii: xml ou legacy
! Variable pour le temps
integer time(maxtim) ! Temps encodes (yyyymmddhh)
integer year(maxtim) ! Annees (yyyy)
integer month(maxtim) ! Mois (mm)
integer day(maxtim) ! Jours (jj)
integer hour(maxtim) ! Heures (hh)
#ifndef NETCDF
integer tstamp(maxtim) ! Liste des "dateo"
#else
real*8 tstamp(maxtim) ! Liste des temps
#endif
character*3 mnth(12) ! Mois de l'annee en abrege
integer ntimes ! Nombre de pas de temps
integer n ! Compteur pour les pas de temps
integer tstep ! Pas de temps (en heures)
integer npas ! Variable FST npas
#ifndef NETCDF
integer dateo ! Variable FST dateo
#else
real*8 dateo
#endif
integer hi ! Heure initiale de la prevision
integer hf ! Heure finale de la prevision
integer ip2 ! Variable FST ip2
integer deet ! Variable FST deet
! Variable pour les trajectoires
real trk_lat(maxctr,maxtrk) ! Latitude
real trk_lon(maxctr,maxtrk) ! Longitude
real trk_lat_f(maxctr,maxtrk) ! Latitude
real trk_lon_f(maxctr,maxtrk) ! Longitude
real trk_vort(maxctr,maxtrk) ! Tourbillon
real trk_press(maxctr,maxtrk) ! Pression
real trk_u(maxctr,maxtrk) ! Estimee pour la vit. de deplacement
real trk_v(maxctr,maxtrk) ! Estimee pour la vit. de deplacement
real trk_circ(maxctr,maxtrk) ! Circulation ou Classe
real trk_b(maxctr,maxtrk) ! Baroclinicite de bas niveaux
real trk_vtl(maxctr,maxtrk) ! Structure de bas niveaux
real trk_vtu(maxctr,maxtrk) ! Structure de hauts niveaux
real trk_maxqr(maxctr,maxtrk) ! Max de tourbillon a 850 hPa
real trk_maxdz(maxctr,maxtrk) ! Max d'epaisseur 250-850 hPa
real trk_maxuv(maxctr,maxtrk) ! Max de vent en surface
real, dimension (:,:) , allocatable :: pmatch ! Probabilite d'association
integer trk_time(maxctr,maxtrk) ! Indice du temps (n)
integer trk_nctr(maxctr) ! Nombre de centres d'une trajectoire
logical trk_closed(maxctr,maxtrk) ! Systeme ferme ou ouvert
logical trk_used(maxtrk) ! Trajectoire active ou terminee
logical trk_warm(maxtrk) ! Trajectoire initialise par un WarmStart
real lat_est ! Position (Latitude) estimee
real lon_est ! Position (Longitude) estimeee
real p_est ! Pression centrale estimee
real v_est ! Tourbillon central estimeee
real prob ! Probabilite d'association trk-ctr
integer it ! Compteur pour les trajectoires
integer itt ! Compteur pour les trajectoires
integer n_track ! Nombre de trajectoires actives
integer n_track_tot ! Nombre de trajectoires totales
integer imatch ! Compteur pour les associations trk-ctr
integer nn ! Nombre de centres d'une trajectoire
integer n_matches ! Nombre d'associations probables trk-ctr
! Variables pour le filtrage (filtre de Cressman)
real srad ! Rayon (m) de filtrage pour le tourbillon
real srad2 ! Rayon (m) de filtrage pour PNM et GZ
real sraduv ! Rayon (m) de filtrage pour les vents
! Variables reliees aux donnees modeles 2D (ni*nj) et 3D(ni*nj*nk)
real, dimension (:,:,:), allocatable :: h3d ! Hauteur multi-niveau
real, dimension (:,:) , allocatable :: h ! Hauteur ou Epaisseur
real, dimension (:,:) , allocatable :: hnf ! Hauteur non filtre
real, dimension (:,:) , allocatable :: hl ! Epaisseur bas niveau
real, dimension (:,:) , allocatable :: u ! Estimation de la vit.
real, dimension (:,:) , allocatable :: v ! de deplacement en x et y
real, dimension (:,:) , allocatable :: press ! PNM
real, dimension (:,:) , allocatable :: pressnf ! PNM non filtre
real, dimension (:,:) , allocatable :: vort ! Tourbillon
real, dimension (:,:) , allocatable :: vct ! Tourbillon avec masque
real, dimension (:,:) , allocatable :: vcrit ! Valeur seuil pour vort
real, dimension (:,:) , allocatable :: orog ! Orographie
real, dimension (:,:) , allocatable :: high_terrain ! Variable pour vcrit
real, dimension (:,:) , allocatable :: windspd ! Vitesse du vent
real, dimension (:,:) , allocatable :: utl ! Vent thermique
real, dimension (:,:) , allocatable :: vtl ! Vent thermique
real, dimension (:,:) , allocatable :: landsea ! Masque Terre-Mer
real, dimension (:,:) , allocatable :: lat ! Latitude
real, dimension (:,:) , allocatable :: lon ! Longitude
real, dimension (:,:) , allocatable :: delta_x ! Distance entre les pts
real, dimension (:,:) , allocatable :: delta_y ! de grille en x et y
real, dimension (:,:) , allocatable :: work1 ! Vecteur de travail
real, dimension (:,:) , allocatable :: work2 ! Vecteur de travail
real, dimension (:,:) , allocatable :: work3 ! Vecteur de travail
real*8, dimension (:,:) , allocatable :: press8 ! PNM
real*8,dimension (:,:) , allocatable :: pressnf8 ! PNM non filtre
real*8, dimension (:,:) , allocatable :: vort8 ! Tourbillon
real*8, dimension (:,:) , allocatable :: h8 ! Hauteur ou Epaisseur
real*8, dimension (:,:) , allocatable :: work8 ! Vecteur de travail
! Variables pour la grille etendue (niext*njext)
real, dimension (:,:) , allocatable :: bufext ! Champ etendu
real, dimension (:,:) , allocatable :: lonext ! Latitude
real, dimension (:,:) , allocatable :: latext ! Longitude
real, dimension (:,:) , allocatable :: delta_xext ! Distance entre les pts
real, dimension (:,:) , allocatable :: delta_yext ! de grille en x et y
real*8, dimension (:,:) , allocatable :: bufext8 ! Champ etendu
! Variables pour le traitement des centres de pression
real, dimension (:) , allocatable :: p_lat ! Latitude
real, dimension (:) , allocatable :: p_lon ! Longitude
real, dimension (:) , allocatable :: p_lat_f ! Latitude
real, dimension (:) , allocatable :: p_lon_f ! Longitude
real, dimension (:) , allocatable :: p_xf_f ! Point de grille en x
real, dimension (:) , allocatable :: p_yf_f ! Point de grille en y
! Variables pour le traitement des centres de tourbillon
real, dimension (:) , allocatable :: v_lat ! Latitude
real, dimension (:) , allocatable :: v_lon ! Longitude
real, dimension (:) , allocatable :: v_xf ! Point de grille en x
real, dimension (:) , allocatable :: v_yf ! Point de grille en y
! Variables pour le traitement des centres de tourbillon
real, dimension (:) , allocatable :: z_lat ! Latitude
real, dimension (:) , allocatable :: z_lon ! Longitude
! Variables pour le traitement des centres suivis
real, dimension (:) , allocatable :: c_vort ! Tourbillon
real, dimension (:) , allocatable :: c_press ! PNM
real, dimension (:) , allocatable :: c_lat ! Latitude
real, dimension (:) , allocatable :: c_lon ! Longitude
real, dimension (:) , allocatable :: c_lat_f ! Latitude
real, dimension (:) , allocatable :: c_lon_f ! Longitude
real, dimension (:) , allocatable :: c_xf ! Point de grille en x
real, dimension (:) , allocatable :: c_yf ! Point de grille en y
real, dimension (:) , allocatable :: c_circ ! Circulation ou Classe
real, dimension (:) , allocatable :: c_u ! Estimation de la vit.
real, dimension (:) , allocatable :: c_v ! de deplacement en x et y
real, dimension (:) , allocatable :: c_b ! Baroclinicite de bas niveau
real, dimension (:) , allocatable :: c_vtl ! Structure de bas niveau
real, dimension (:) , allocatable :: c_vtu ! Structure de haut niveau
real, dimension (:) , allocatable :: c_maxqr ! Maximum de tourbillon
real, dimension (:) , allocatable :: c_maxdz ! Maximum d'epaiseur
real, dimension (:) , allocatable :: c_maxuv ! Vent de surface maximum
real, dimension (:) , allocatable :: c_landsea ! Vent de surface maximum
real, dimension (:) , allocatable :: wc_lat ! Latitude - WarmStart
real, dimension (:) , allocatable :: wc_lon ! Longitude - WarmStart
logical, dimension (:) , allocatable :: c_warm ! Centre actif au temp n = 1
logical, dimension (:) , allocatable :: c_closed ! Systeme ouvert ou ferme
logical, dimension (:) , allocatable :: c_used ! Centre utilise ou non
integer, dimension (:) , allocatable :: icm ! Centres avec prob > 0
integer, dimension (:) , allocatable :: itm ! Tracks associes a icm
integer, dimension (:) ,allocatable :: nivh ! Niv. de geopot. interpole
integer nivhm(K_NIVMAX) ! Nombre de niveau de geopotentiel en entree
integer nc ! Nombre de centres max
integer iv ! Compteur des centres de tourbillon
integer ip_f ! Compteur des centres de pression filtree
integer ip ! Compteur des centres de pression brute
integer iz ! Compteur des centres d'epaisseur
integer ic ! Compteur pour les centres suivis
integer n_new ! Nombre de nouveaux centres detectes
integer n_warm ! Nombre de centre dans le fichier ./warmstart.txt
integer ihem ! Facteur pour la conversion du tourbillon
integer nkh ! Nombre de niveaux de GZ interpolle
integer nkhm ! Nombre de niveaux de GZ en entree
real vort_t ! Valeur interpollee du tourbillon
real dist ! Distance (m)
real dummy ! Valeur sans usage
real f_lat ! Latitude du centre le plus pres
real f_lon ! Longitude du centre le plus pres
logical found ! True = un extreme ce trouve en ce point i,j
logical start ! True = on demarre une trajectoire
! Variables associees aux fichers d'entrees/sorties
character*1024 infile ! Fichier d'entree
character*1024 infile2 ! Fichier d'entree orographie
character*1024 infile3 ! Fichier d'entree masque terre-mer
character*1024 outfile ! Base pour les fichier de sorties
character*1024 ofname1 ! Fichier de sortie #1
character*1024 ofname2 ! Fichier de sortie #2
character*1024 configfile ! Fichier de configuration json
integer iun01 ! Unite pour le fichier d'entree
integer iun02 ! Unite pour le fichier d'entree orographie
integer iun03 ! Unite pour le fichier d'entree masque terre-mer
integer iun51 ! Unite pour le fichier de sortie #1
integer iun52 ! Unite pour le fichier de sortie #2
integer iun60 ! Unite pour le fichier ./warmstart.txt
integer real_len ! Fonction pour la calcul de la longueur d'un string
#ifndef NETCDF
integer fnom ! Fonction ARMNLIB
integer fstouv ! Fonction ARMNLIB
integer fstopc ! Fonction ARMNLIB
integer fstfrm ! Fonction ARMNLIB
integer fclos ! Fonction ARMNLIB
#else
integer wstart(4), wcount(4)
integer ilatid, ilonid, itimeid, latdimid, londimid, timedimid
integer lonextdimid, latextdimid
integer latid, lonid, timeid
#endif
integer nrecs(3) ! Valeur retournee par FSTOUV
integer ier ! Valeur retournee par FSTLIR
logical warmstart ! Vrai si un fichier ./warmstart.txt existe
! Variables pour la manipulation des grilles
real xpf ! Point de grille interpole en x
real ypf ! Point de grille interpole en y
#ifndef NETCDF
character*2 getvar ! Variable pour l'obtention des parametres de grille
character*1 typvar ! Type des variable RPN
character*1 grtyp ! Type de grille
#else
character*256 getvar ! Variable pour l'obtention des parametres de grille
character*256 attname
character*256 attvalue
integer rhdims(4)
integer orogid, orogfid
character*(NF_MAX_NAME) latname, lonname
#endif
integer global ! 0 = grille LAM, 1 = grille globale ou ni /= i
! 2 = grille globale ou ni=i (repetition)
integer istart ! Point de grille de depart en i
integer jstart ! Point de grille terminal en j
integer iend ! Point de grille de depart en i
integer jend ! Point de grille terminal en j
#ifndef NETCDF
integer gdid ! Descripteur de grille pour EZSCINT
integer gdidext ! Descripteur de la grille etendue pour EZSCINT
integer ezqkdef ! Fonction EZSCINT
#endif
integer gdll ! Fonction EZSCINT
integer gdxysval ! Fonction EZSCINT
integer gdllsval ! Fonction EZSCINT
integer gdxyfll ! Fonction EZSCINT
integer gdllfxy ! Fonction EZSCINT
integer i ! Compteur pour les points de grille en x
integer j ! Compteur pour les points de grille en y
integer ni ! Dimension de la grille en x
integer nj ! Dimension de la grille en y
integer ip1 ! Variable FST ip1
integer ip3 ! Variable FST ip3
integer glbkey ! Descripteur pour la grille
integer glbkeygz ! Descripteur pour les niveaux de GZ
integer get_dimgrille ! Fonction pour les dimensions de la grille
#ifndef NETCDF
integer fstprm_dio ! Fonction pour les parametres de grille
#endif
integer ig1 ! Descripteur FST de la grille d'origine
integer ig2 ! " "
integer ig3 ! " "
integer ig4 ! " "
integer ig1ext ! Descripteur FST de la grille entedue
integer ig2ext ! " "
integer ig3ext ! " "
integer ig4ext ! " "
integer niext ! Dimension de la grille etendue en x
integer njext ! Dimension de la grille etendue en y
! Variable de configuration des entrees
type :: config_type
character(len=:),allocatable :: generic_name
character(len=:),allocatable :: name
character(len=:),allocatable :: units_destination
real :: level
real :: scale
real :: offset
end type config_type
type(config_type), dimension(13) :: config
integer nconfig, curn
character*20 str
character*30 tmpstr
character*256 strerr
!**********************************************************************
! Traitement des cles d'appel
!**********************************************************************
call get_options(SILENT, IP3, GETVAR,
$ INFILE, INFILE2, INFILE3, OUTFILE, CONFIGFILE,
$ NPT, DIST_CRIT, LEVEL, FVORT, VC, ANTICYCL, VORTCNT,
$ ANALYSE, TROPIC, NOFST, FPRES, TCRAD, HCRIT, TDWIND,
$ PHASE, BCRIT, LATMIN, LATMAX, LONMIN, LONMAX, RELOCP,
$ ONLYWARM, OUTPUT_TYPE)
!**********************************************************************
! Initialisations
!**********************************************************************
data mnth/'JAN','FEB','MAR','APR','MAY','JUN',
$ 'JUL','AUG','SEP','OCT','NOV','DEC'/
iun01 = 1
iun02 = 2
iun03 = 3
iun51 = 51
iun52 = 52
iun60 = 60
#ifdef NETCDF
! Fichier de configuration JSON
call read_config(configfile, config, nconfig)
#endif
#ifndef NETCDF
ofname1 = outfile(1:real_len(outfile))//'.fst' ! fichier sortie RPN
#else
ofname1 = outfile(1:real_len(outfile))//'.nc' ! fichier sortie NetCDF
#endif
if ( trim(output_type) == 'xml' ) then
ofname2 = outfile(1:real_len(outfile))//'.xml' ! fichier sortie XML
else
ofname2 = outfile(1:real_len(outfile))//'.txt' ! fichier sortie TXT
endif
#ifndef NETCDF
! Ouverture fichier input RPN
ier = fnom(iun01, infile(1:real_len(infile)), 'RND', 0)
nrecs(1) = fstouv(iun01, 'RND')
if (nofst.ne.0) then
! Ouverture fichier output RPN
ier = fnom(iun51, ofname1, 'RND', 0)
nrecs(2) = fstouv(iun51, 'RND')
endif
#else
! Ouverture fichier input NETCDF
ier = nf_open(infile(1:real_len(infile)), NF_NOWRITE, iun01)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error opening NetCDF input file: ',infile(1:real_len(infile))
call handle_err(ier,strerr)
endif
! Ouverture fichier input NETCDF orographie
ier = nf_open(infile2(1:real_len(infile2)), NF_NOWRITE, iun02)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error opening NetCDF input file 2: ',infile2(1:real_len(infile2))
call handle_err(ier,strerr)
endif
! Ouverture fichier input NETCDF masque terre-mer
ier = nf_open(infile3(1:real_len(infile3)), NF_NOWRITE, iun03)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error opening NetCDF input file 3: ',infile3(1:real_len(infile3))
call handle_err(ier,strerr)
endif
if (nofst.ne.0) then
! Ouverture fichier output NetCDF
ier = nf_create(ofname1, OR(NF_CLASSIC_MODEL, NF_NETCDF4), iun51)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error opening NetCDF output file: ',ofname1(1:real_len(ofname1))
call handle_err(ier,strerr)
endif
endif
#endif
! Ouverture du fichier des trajectoires (txt)
open(unit=iun52, file=ofname2, status='unknown', form='FORMATTED')
if ( trim(output_type) == 'xml' ) then
write(iun52,'(a)') "<tracks>"
endif
#ifndef NETCDF
! FSTD: On ne tolere pas les erreurs a partir de warning et plus
ier = fstopc('TOLRNC','WARNIN',.false.)
#endif
!**********************************************************************
! Parametres des grilles
!**********************************************************************
! Dimension de la grille en entree
glbkey = get_dimgrille(ni, nj, iun01, getvar)
#ifndef NETCDF
! Parametres pour l'ecriture en fichier standard RPN
ier=fstprm_dio(dateo,deet,typvar,grtyp,ig1,ig2,ig3,ig4,glbkey)
! Parametres pour EZSCINT
gdid = ezqkdef(ni,nj,grtyp,ig1,ig2,ig3,ig4,iun01)
#else
! Parametres meta-donnees globales fichier NETCDF
if (nofst.ne.0) then
attname = "Conventions"
attvalue = "CF-1.6"
ier = nf_put_att_text(iun51, NF_GLOBAL, attname, real_len(attvalue),
$ attvalue(1:real_len(attvalue)))
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error writing NetCDF attribute to output file ',
$ ofname1(1:real_len(ofname1)),' : ',attname,'=',attvalue(1:real_len(attvalue))
call handle_err(ier,strerr)
endif
attname = "title"
attvalue = "Sinclair Tracking diagnostic file"
ier = nf_put_att_text(iun51, NF_GLOBAL, attname, real_len(attvalue),
$ attvalue(1:real_len(attvalue)))
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error writing NetCDF attribute to output file ',
$ ofname1(1:real_len(ofname1)), ' : ',attname,'=',attvalue(1:real_len(attvalue))
call handle_err(ier,strerr)
endif
endif
#endif
! Nombre de niveau de geopotentiel
if (tropic .ge. 0 .or. phase .ge. 0) then
nkh = 13
endif
!**********************************************************************
! Allocation de memoire dynamique
!**********************************************************************
! Grille d'origine
allocate ( h(ni,nj), hnf(ni,nj), u(ni,nj), v(ni,nj) )
allocate ( press(ni,nj), pressnf(ni,nj) )
allocate ( vort(ni,nj), vct(ni,nj), vcrit(ni,nj) )
allocate ( windspd(ni,nj))
allocate ( orog(ni,nj), high_terrain(ni,nj) )
allocate ( lat(ni,nj), lon(ni,nj) )
allocate ( delta_x(ni,nj), delta_y(ni,nj) )
allocate ( work1(ni,nj), work2(ni,nj), work3(ni,nj) )
allocate ( press8(ni,nj), vort8(ni,nj), h8(ni,nj), work8(ni,nj) )
allocate ( pressnf8(ni,nj))
! Centres
nc = nint(ni*nj/4.0)
allocate ( p_lat(nc), p_lon(nc) )
allocate ( p_lat_f(nc), p_lon_f(nc) )
allocate ( p_xf_f(nc), p_yf_f(nc) )
allocate ( v_lat(nc), v_lon(nc) )
allocate ( v_xf(nc), v_yf(nc) )
if (tropic .ge. 0 .or. phase .ge. 0) then
allocate ( nivh(nkh) )
allocate ( h3d(ni,nj,nkh) )
allocate ( c_b(nc) )
allocate ( c_maxuv(nc) )
allocate ( hl(ni,nj), utl(ni,nj), vtl(ni,nj))
allocate ( z_lat(nc), z_lon(nc) )
allocate ( c_maxqr(nc), c_maxdz(nc) )
allocate ( c_vtl(nc), c_vtu(nc) )
endif
if (tropic .ge. 0) then
allocate ( landsea(ni,nj) )
allocate ( c_landsea(nc))
endif
allocate ( c_vort(nc), c_press(nc), c_lat(nc), c_lon(nc) )
allocate ( c_lat_f(nc), c_lon_f(nc), c_warm(nc) )
allocate ( c_xf(nc), c_yf(nc) )
allocate ( c_circ(nc), c_u(nc), c_v(nc), c_closed(nc) )
allocate ( c_used(nc), icm(nc), itm(nc) )
! Probabilites d'association
allocate ( pmatch(nc,maxtrk) )
!**********************************************************************
! Latitude-Longitude et distance entre les points de grille
!**********************************************************************
#ifndef NETCDF
ier = gdll(gdid, LAT, LON)
#else
ier = gdll(LAT, LON, LATNAME, LONNAME, iun01, glbkey, ni, nj)
#endif
#ifndef NETCDF
if (nofst.ne.0) then
call ecrit_fst(iun51,lat,'LA','LATITUDE',0,0,
$ 0,dateo,ni,nj,ig1,ig2,ig3,ig4,ip3,grtyp,typvar,deet)
call ecrit_fst(iun51,lon,'LO','LONGITUD',0,0,
$ 0,dateo,ni,nj,ig1,ig2,ig3,ig4,ip3,grtyp,typvar,deet)
endif
#endif
c Grille globale ou aire limitee ?
call global_or_lam(GLOBAL, lon, ni, nj)
c On evalue la distance entre les points de grille
call calc_dxdy_gr(DELTA_X, DELTA_Y, lat, lon, ni, nj)
!**********************************************************************
! Grille etendue pour le filtrage
!**********************************************************************
#ifndef NETCDF
if ( global .eq. 0 ) then ! LAM
call extend_lam_grid(NIEXT, NJEXT, ig1ext, ig2ext, ig3ext,
$ ig4ext, glbkey, ni, nj)
gdidext = ezqkdef(niext, njext, grtyp, ig1ext, ig2ext, ig3ext,
$ ig4ext, iun01)
else
call extend_glb_grid(NIEXT, glbkey, ni)
njext = nj
ig1ext = ig1
ig2ext = ig2
ig3ext = ig3
ig4ext = ig4
endif
#else
if ( global .eq. 0 ) then ! LAM
call extend_lam_grid(NIEXT, NJEXT, lat, lon, ni)
else
call extend_glb_grid(NIEXT, lon, ni, nj)
njext = nj
ig1ext = ig1
ig2ext = ig2
ig3ext = ig3
ig4ext = ig4
endif
#endif
allocate ( bufext(niext,njext) )
allocate ( latext(niext,njext), lonext(niext,njext) )
allocate ( delta_xext(niext,njext), delta_yext(niext,njext) )
allocate ( bufext8(niext,njext) )
if ( global .eq. 0 ) then ! LAM
#ifndef NETCDF
ier = gdll(gdidext, LATEXT, LONEXT)
#endif
call calc_dxdy_gr(DELTA_XEXT, DELTA_YEXT, latext, lonext,
$ niext, njext)
else
call extend_glb(DELTA_XEXT, delta_x, global, ni, nj, niext)
call extend_glb(DELTA_YEXT, delta_y, global, ni, nj, niext)
endif
#ifdef NETCDF
ier = nf_def_dim(iun51, 'longitudeext', niext, lonextdimid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error defining NetCDF dimension in output file ',
$ ofname1(1:real_len(ofname1)), ' : longitudeext'
call handle_err(ier,strerr)
endif
ier = nf_def_dim(iun51, 'latitudeext', njext, latextdimid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error defining NetCDF dimension in output file ',
$ ofname1(1:real_len(ofname1)), ' : latitudeext'
call handle_err(ier,strerr)
endif
#endif
!**********************************************************************
! Lecture des parametres lies au temps
!**********************************************************************
#ifndef NETCDF
call get_times(YEAR, MONTH, DAY, HOUR, TSTAMP, NTIMES, TSTEP,
$ HI, HF, iun01, 'GZ', level, maxtim)
#else
call get_times(YEAR, MONTH, DAY, HOUR, TSTAMP, NTIMES, TSTEP,
$ iun01, 'time', maxtim)
#endif
!**********************************************************************
! On regle certaines options de l'algorithme
!**********************************************************************
call set_track_options(LEVEL, ISIGN, NPT, DIST_CRIT, VORT_C,
$ DO_CIRC, FVORT, FPRES, VC, TCRAD, HCRIT, TDWIND, BCRIT,
$ vortcnt, anticycl, tropic, tstep, phase)
if ( tropic .ge. 0 .or. phase .ge. 0) then
#ifndef NETCDF
! On trouve le niveau qui represente le vent a 10m
call find_uvsurflev(SLEVEL,iun01,ni,nj)
#endif
#ifndef NETCDF
! On cherche les niveaux de pression de GZ du fichier en entree
glbkeygz = get_dimgrille(ni, nj, iun01, 'GZ')
call get_niv(nivhm, nkhm, glbkeygz, iun01)
#else
call get_config(curn, config, "height", nconfig)
if (curn .gt. 0) then
tmpstr = config(curn)%name
else
tmpstr = "zg"
endif
call get_niv(nivhm, nkhm, tmpstr, iun01)
#endif
print*
print*,"Niveaux (hPa) de la hauteur du geopotentiel en entree"
print*,nkhm, nivhm(1:nkhm)
! On definie les niveaux de la hauteur du geopotentiel pour l'interpolation verticale
nivh(1)=300 ! hPa
nivh(2)=350
nivh(3)=400
nivh(4)=450
nivh(5)=500
nivh(6)=550
nivh(7)=600
nivh(8)=650
nivh(9)=700
nivh(10)=750
nivh(11)=800
nivh(12)=850
nivh(13)=900
endif
!**********************************************************************
! Ecriture header + dimensions fichier NetCDF si ecriture fichier de sortie
!**********************************************************************
#ifdef NETCDF
if (nofst.ne.0) then
! Dimensions
ier = nf_def_dim(iun51, 'longitude', ni, londimid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error defining NetCDF dimension in output file ',
$ ofname1(1:real_len(ofname1)), ' : longitude'
call handle_err(ier,strerr)
endif
ier = nf_def_dim(iun51, 'latitude', nj, latdimid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error defining NetCDF dimension in output file ',
$ ofname1(1:real_len(ofname1)), ' : latitude'
call handle_err(ier,strerr)
endif
ier = nf_def_dim(iun51, 'time', NF_UNLIMITED, timedimid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error defining NetCDF dimension in output file ',
$ ofname1(1:real_len(ofname1)), ' : time'
call handle_err(ier,strerr)
endif
! Define dimension variables
rhdims(1) = timedimid
ier = nf_def_var(iun51, 'time', nf_double, 1, rhdims, timeid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error defining NetCDF variable in output file ',
$ ofname1(1:real_len(ofname1)), ' : time'
call handle_err(ier,strerr)
endif
rhdims(1) = londimid
rhdims(2) = latdimid
ier = nf_def_var(iun51, 'longitude', nf_float, 2, rhdims, lonid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error defining NetCDF variable in output file ',
$ ofname1(1:real_len(ofname1)), ' : longitude'
call handle_err(ier,strerr)
endif
rhdims(1) = londimid
rhdims(2) = latdimid
ier = nf_def_var(iun51, 'latitude', nf_float, 2, rhdims, latid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error defining NetCDF variable in output file ',
$ ofname1(1:real_len(ofname1)), ' : latitude'
call handle_err(ier,strerr)
endif
! Copy CF attributes from input file
ier = nf_inq_varid(iun01, 'time', itimeid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error getting NetCDF variable ID in input file ',
$ infile(1:real_len(infile)),' : time'
call handle_err(ier,strerr)
endif
ier = nf_copy_att(iun01, itimeid, 'standard_name', iun51, timeid)
if (ier .ne. nf_noerr) then
print*,"Warning: No standard_name attribute for time for variable time."
endif
ier = nf_copy_att(iun01, itimeid, 'long_name', iun51, timeid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error copying NetCDF attribute ID from input file ',infile(1:real_len(infile)),
$ ' to output file ', ofname1(1:real_len(ofname1)), ' for variable time: long_name'
call handle_err(ier,strerr)
endif
ier = nf_copy_att(iun01, itimeid, 'units', iun51, timeid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error copying NetCDF attribute ID from input file ',infile(1:real_len(infile)),
$ ' to output file ',ofname1(1:real_len(ofname1)),' for variable time: units'
call handle_err(ier,strerr)
endif
ier = nf_copy_att(iun01, itimeid, 'calendar', iun51, timeid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error copying NetCDF attribute ID from input file ',infile(1:real_len(infile)),
$ ' to output file ',ofname1(1:real_len(ofname1)),' for variable time: calendar'
call handle_err(ier,strerr)
endif
ier = nf_inq_varid(iun01, latname, ilatid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error getting NetCDF variable ID in input file ',
$ infile(1:real_len(infile)),' : latitude'
call handle_err(ier,strerr)
endif
ier = nf_copy_att(iun01, ilatid, 'standard_name', iun51, latid)
if (ier .ne. nf_noerr) then
print*,"Warning: No standard_name attribute for latitude."
endif
ier = nf_copy_att(iun01, ilatid, 'long_name', iun51, latid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error copying NetCDF attribute ID from input file ',
$ infile(1:real_len(infile)),' to output file ',
$ ofname1(1:real_len(ofname1)), ' for variable latitude: long_name'
call handle_err(ier,strerr)
endif
ier = nf_copy_att(iun01, ilatid, 'units', iun51, latid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error copying NetCDF attribute ID from input file ',
$ infile(1:real_len(infile)),' to output file ',
$ ofname1(1:real_len(ofname1)), ' for variable latitude: units'
call handle_err(ier,strerr)
endif
ier = nf_inq_varid(iun01, lonname, ilonid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error getting NetCDF variable ID in input file ',
$ infile(1:real_len(infile)),' : longitude'
call handle_err(ier,strerr)
endif
ier = nf_copy_att(iun01, ilonid, 'standard_name', iun51, lonid)
if (ier .ne. nf_noerr) then
print*,"Warning: No standard_name attribute for longitude."
endif
ier = nf_copy_att(iun01, ilonid, 'long_name', iun51, lonid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error copying NetCDF attribute ID from input file ',
$ infile(1:real_len(infile)),' to output file ',
$ ofname1(1:real_len(ofname1)), ' for variable longitude: long_name'
call handle_err(ier,strerr)
endif
ier = nf_copy_att(iun01, ilonid, 'units', iun51, lonid)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error copying NetCDF attribute ID from input file ',
$ infile(1:real_len(infile)),' to output file ',
$ ofname1(1:real_len(ofname1)), ' for variable longitude: units'
call handle_err(ier,strerr)
endif
! End define mode
ier = nf_enddef(iun51)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error closing NetCDF definition mode for output file ',
$ ofname1(1:real_len(ofname1))
call handle_err(ier,strerr)
endif
! Write dimension variables
wstart(1) = 1
wstart(2) = 1
wcount(1) = ntimes
ier = nf_put_vara_double(iun51, timeid, wstart, wcount, tstamp)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error cannot write time variable to output file ',
$ ofname1(1:real_len(ofname1))
call handle_err(ier,strerr)
endif
wcount(1) = ni
wcount(2) = nj
ier = nf_put_vara_real(iun51, lonid, wstart, wcount, lon)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error cannot write longitude variable to output file ',
$ ofname1(1:real_len(ofname1))
call handle_err(ier,strerr)
endif
wcount(1) = ni
wcount(2) = nj
ier = nf_put_vara_real(iun51, latid, wstart, wcount, lat)
if (ier .ne. nf_noerr) then
write(strerr,*) 'Error cannot write latitude variable to output file ',
$ ofname1(1:real_len(ofname1))
call handle_err(ier,strerr)
endif
endif
call ecrit_nc(iun51, delta_x, 'delta_x', 1, 1,
$ londimid, latdimid, timedimid, ni, nj)
call ecrit_nc(iun51, delta_y, 'delta_y', 1, 1,
$ londimid, latdimid, timedimid, ni, nj)
#endif
!**********************************************************************
! Lecture de la topographie et calcul de VCRIT
!**********************************************************************
#ifndef NETCDF
if (tropic .ge. 0) then
call read_landsea(LANDSEA,iun01,ni,nj)
endif
call read_orog(OROG,iun01,ni,nj)
#else
if (tropic .ge. 0) then
call read_landsea(LANDSEA,iun03,'lsm',ni,nj)
endif
call read_orog(OROG,iun02,'z',ni,nj)
#endif
if (nofst.ne.0) then
#ifndef NETCDF
call ecrit_fst(iun51,orog,'MX','TOPO_NF ',0,0,
$ 0,dateo,ni,nj,ig1,ig2,ig3,ig4,ip3,grtyp,typvar,deet)
#else
call ecrit_nc(iun51, orog, 'z', 1, 1, londimid, latdimid,
$ timedimid, ni, nj)
#endif
endif
! smooth orography to 400. km
work8(:,:) = real(orog(:,:),kind=8)
call smooth_c(WORK8, ni, nj, 400.E3, work8, delta_x, delta_y)
orog(:,:) = real(work8(:,:))
if (nofst.ne.0) then
#ifndef NETCDF
call ecrit_fst(iun51,orog,'MX','TOPO ',0,0,
$ 0,dateo,ni,nj,ig1,ig2,ig3,ig4,ip3,grtyp,typvar,deet)
#else
call ecrit_nc(iun51, orog, 'zf', 1, 1, londimid, latdimid,
$ timedimid, ni, nj)
endif
#endif
! Pour le suivi pres de la surface, on ajuste la valeur critique
! du tourbillon (vcrit) en fonction de la topographie
if (vortcnt.eq.0 .and. level.ge.850) then
call calc_vcrit(VCRIT, HIGH_TERRAIN, orog, vc, fvort, ni, nj)
else
VCRIT(:,:) = vc / 1.E5
HIGH_TERRAIN(:,:) = 0.0
endif
!**********************************************************************
! WarmStart
!**********************************************************************
! On regarde si un fichier ./warmstart.txt existe
inquire(file="./warmstart.txt", exist=warmstart)
if (warmstart) then
print*
print*, "mode WarmStart ACTIVE"
open(unit=iun60, file="./warmstart.txt",action="read",
$ status="old")
n_warm = 0
do
read (iun60,*, end=10)
n_warm = n_warm + 1
end do
10 close (iun60)
if (n_warm .ge. 1) then
allocate(wc_lat(n_warm))