Philippe 23/02/2016: lfi2cdf: modif for JPHEXT/=1
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / EXTRACTDIA / obs2mesonh.f90
1 !     ###################
2       PROGRAM  OBS2MESONH
3 !     ###################
4 !
5 !!****  *OBS2MESONH* -  Interpolation d une liste de valeurs observées
6 !                       sur la grille Mesonh
7 !                       en entrée un fichier ASCII au format [ ] pour optionnel
8 !                                 [YYYMMJJHHMMSS]
9 !                                 lon lat [alt] valeur_obs  
10 !                              ou lat lon [alt] valeur_obs 
11 !                       en sortie un fichier diachronique
12 !! 
13 !!
14 !!    PURPOSE
15 !!    -------
16
17 !
18 !!**  METHOD
19 !!    ------
20 !!      
21 !     Lecture en entree:
22 !  d un fichier ascii contenant les localisations (lon,lat alt,valeur) a traiter
23 !  d un fichier modele pour recuperer la grille XYZ
24 !     Pour chaque obs lue, 1) recherche du point de grille xy la contenant,
25 !                          2) recherche sur la verticale du niveau K la contenant
26 !     Si plusieurs obs sont contenues dans un meme point de grille, calcul de la moyenne des ces obs
27 !     Pour certaines variables (unite dBz, nom de champ_WVBT ou _IRBT), 
28 !   passage a des unites plus pertinentes pour effectuer les moyennes
29 !   et retour aux unites d origine avant ecriture (voir les 2 routines
30 ! symetriques To_computing_units et From_computing_units)
31 !     Mise a XSPVAL des points de grille ne contenant pas d'obs
32 !      
33 !     Ecriture en sortie:
34 !  d un fichier  diachronique ( utiliser LSPOT=T dans diaprog pour visualiser 
35 !       toutes les points grilles non XSPVAL)
36 !
37 !!
38 !!    EXTERNAL
39 !!    --------
40 !!          CREATLINK : à l'ouverture du fichier, HYFLAGFILE='OPE',
41 !!                      création d'un lien dans le directory local
42 !!                      si le fichier existe sous $DIROBS
43 !!          READVAR   : lecture d unchamp du fichier diachronique
44 !!          WRITEVAR : ecriture format lon lat alt val
45 !!          SM_XYHAT  : création de la grille des Obs
46 !!          TO_COMPUTING_UNITS: passage unites vers unites plus pertinentes 
47 !!                              pour effectuer des calculs
48 !!          FROM_COMPUTING_UNITS: passage inverse avant ecriture
49 !!                               (appele par writevar)
50 !!
51 !!    REFERENCE
52 !!    ---------
53 !!
54 !!    AUTHORS
55 !!    -------
56 !!    N. Asencio and J. Stein
57 !!
58 !!
59 !!    MODIFICATIONS
60 !!    -------------
61 !!      Original    04/11/2003
62 !!      Fev 2005: ajout de champs diachroniques ALT_champ N_champ
63 !!                changement de grille pour le vent (zonal,meridien->
64 !!                grille Mesonh)
65 !!      04/05/2005 add a control for the min and max of the field before
66 !                and after interpolation(s)
67 !                  observations outside the mesonh domain are rejected
68 !       19/09/2005 G.Jaubert CNRM  
69 !                  1) Nom du fichier .lfi en output demande 
70 !                  2) l'enregistrement peut ne pas contenir alt si champ 2D
71 !                  3) si le fichier de donnees commence par une date,
72 !                     reinitialisation de la date dans le lfi de sortie
73 !-------------------------------------------------------------------------------
74 !
75 !*       0.    DECLARATIONS
76 !              ------------
77 #ifdef NAGf95
78 USE F90_UNIX  ! for FLUSH
79 USE F90_UNIX_PROC  ! for SYSTEM
80 #endif      
81 ! modules MesoNH
82 USE MODD_CST
83 USE MODD_PARAMETERS, ONLY:JPHEXT,JPVEXT,XUNDEF
84 USE MODD_DIM1, ONLY:NIMAX,NJMAX,NKMAX, NIINF, NISUP ,NJINF,NJSUP
85 USE MODD_GRID, ONLY: XLON0,XLAT0,XLONORI,XLATORI
86 USE MODD_GRID1, ONLY: XXHAT,XYHAT,XZZ
87 USE MODE_GRIDPROJ  ! subroutine SM_XYHAT 
88 ! modules DIACHRO
89 USE MODN_NCAR,  ONLY: XSPVAL      
90 USE MODD_COORD
91 !                    XVAR(i,j,k,,,), XMASK,XTRAJ ,XDATIME(16,t)
92 !                     et NGRIDIA , NGRIDIAM ( appel interp_grids)
93 USE MODD_ALLOC_FORDIACHRO
94 !                      nverbia, CGROUP
95 USE MODD_RESOLVCAR 
96 USE MODD_READLH !NREADIL,IH,...
97 !
98 !
99 USE MODI_WRITEVAR
100 USE MODI_CREATLINK
101 USE MODI_LOW2UP      
102 USE MODI_WRITEDIR      
103 USE MODI_UV_TO_ZONAL_AND_MERID
104 USE MODI_TO_COMPUTING_UNITS
105 !
106 IMPLICIT NONE                       
107 !
108 !*       0.1   Local variables declaration
109 !
110 !! indices de boucle
111 INTEGER     :: JILOOP,JOBS
112 ! zoom  suivant les 6 dimensions des champs diachro
113 INTEGER     :: ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin
114 INTEGER     :: ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup
115 INTEGER     :: IL   ! indice de positionnement dans yligne
116 INTEGER     :: IIMNH,IJMNH,IKMNH,IKMAX,IGRID
117 !
118 INTEGER     :: IAN, IMOIS, IJOUR, IHEUR, IMINU, ISEC ! date observation
119 REAL        :: XSEC ! heure observation en secondes 
120 !                                   
121 ! stockage
122 REAL, allocatable, dimension(:,:,:)    ::  ZOBSinMNH,ZALTOBS
123 INTEGER , allocatable, dimension(:,:,:) ::  ICPTOBSinMNH    
124 REAL              :: ZXOBS,ZYOBS,ZOBSLONLU,ZOBSLATLU,ZOBSALTLU,ZVALOBS
125 ! pour le passage composantes meridienne,zonale a la grille Mesonh
126 REAL, allocatable, dimension(:,:,:,:,:,:):: ZVENTSAVE
127 REAL, allocatable, dimension(:,:,:)      :: ZWORK3D,ZWORK3D2
128 !
129 REAL :: zmini,zmaxi
130 REAL :: ZVAR1, ZVAR2, ZVAR3
131 LOGICAL :: galtobs
132 INTEGER           :: ILUDIR,iret,ilocverbia,inbval,inbvalrej,IKM1
133 !! **** la taille des variables caracteres contenant les noms
134 !!      de fichiers est obligatoirement de 28 ****
135 !!      pour toutes les routines diachro
136 CHARACTER(LEN=28) :: YFILEGRID,YDUMMYFILE , YFILEOUTNAME
137 CHARACTER(LEN=100):: YFILEOBS
138 CHARACTER(LEN=3)  :: YFLAGREADVAR ,YFLAGWRITE
139 CHARACTER (LEN=10):: YUNITE,YUNITEMAJ
140 CHARACTER (LEN=3) :: YSTOCK,YFILEOUTSUFFIX
141 CHARACTER (LEN=4) :: YLL
142 CHARACTER(LEN=11) :: YLUDIR      !  Name of the dir file
143 CHARACTER(LEN=14) :: YDATEOBS    ! observation  date (YYYYMMDDHHMISS) 
144 CHARACTER(LEN=100) :: YLIGNE    ! 
145 !-------------------------------------------------------------------------------
146 !
147 !*       1.    Init
148 !              ----
149 !
150 ! active(1) ou desactive(0) les prints de controle dans les routines
151 ! READVAR et WRITEVAR
152 ilocverbia=0
153
154 ! dans mesonh Xundef est utilise 
155 ! dans les routines diachro XSPVAL est utilisé
156  XSPVAL=XUNDEF                                    
157 !
158 ! ouverture d un fichier dir ou vont s ecrire les entrees clavier
159 YLUDIR='dirobs2mnh'
160 CALL FMATTR(YLUDIR,YLUDIR,ILUDIR,iret)
161 OPEN(UNIT=ILUDIR,FILE=YLUDIR,FORM='FORMATTED')
162 !
163 NIINF=0
164 NISUP=0
165 NJINF=0
166 NJSUP=0
167 iret=0
168 !
169 !*      2.  Lecture et initialisation des modules Mesonh
170 !          ----------------------------
171
172 PRINT*, '- Name of the diachro file to read the grid '&
173       ,'(without .lfi) ?'
174 READ(5,'(A)',END=99) YFILEGRID
175 YFILEGRID=ADJUSTL(YFILEGRID)      
176 CALL WRITEDIR(ILUDIR,YFILEGRID)
177 !
178 PRINT*, '- Prints : 0= mini 1=debug mode in obs2mesonh'
179 PRINT*, '                   2= print of input values'
180 PRINT*, '                   3=debug mode in diachro routines'
181 PRINT*, '?'
182 READ(5,*)ilocverbia
183 CALL WRITEDIR(ILUDIR,ilocverbia)
184 PRINT*, ' output prints= ',ilocverbia 
185 IF (ilocverbia >2) nverbia=ilocverbia
186 !
187 ! Lecture du  champ ZSBIS  pour  obtenir
188 ! l Initialisation des variables
189 ! des modules (cf USE en debut de programme)
190 !
191 !  indique que le fichier lu doit etre ouvert dans READVAR
192 YFLAGREADVAR='OPE'
193 CALL READVAR('ZSBIS',YFILEGRID,YFLAGREADVAR,ilocverbia,iret)
194 print *, 'READVAR(zsbis), return value= ',iret
195 IF ( iret /= 0 ) THEN 
196   print *,'** Error when reading the grid in the FM diachro file: ',TRIM(YFILEGRID)
197   STOP
198 ENDIF
199 !
200 if (ilocverbia >= 0 ) then
201   print *,' Size of input array(zs)= ',SIZE(XVAR,1),SIZE(XVAR,2),&
202             SIZE(XVAR,3),SIZE(XVAR,4),SIZE(XVAR,5),SIZE(XVAR,6)
203   PRINT*, 'NIINF,NISUP,NJINF,NJSUP', NIINF,NISUP,NJINF,NJSUP
204   PRINT*, 'LATORI,LONORI ', XLATORI,XLONORI
205   zmini= 0.5 * (XXHAT(1)+XXHAT(2))
206   zmaxi= 0.5 * (XYHAT(1)+XYHAT(2))
207   CALL SM_LATLON(XLATORI,XLONORI,zmini,zmaxi,ZVAR1,ZVAR2)
208   PRINT*, 'LATOR,LONOR ', ZVAR1,ZVAR2
209 endif
210 !
211 !
212 PRINT*,'- Name of the output file ? (the default [CR or empty line], "obs" will be added to the input FM file'
213 READ(5,'(A28)',END=88) YFILEOUTNAME
214 YFILEOUTNAME=ADJUSTL(YFILEOUTNAME)
215 CALL WRITEDIR(ILUDIR,YFILEOUTNAME)
216 IF (YFILEOUTNAME(1:2) == 'll' .OR. YFILEOUTNAME(1:2) == 'LL' ) THEN
217       print * ,'** OBS2MESONH: the format of the input file was modified Oct2005'
218       print * ,' the 3rd line is for the name of the output file'
219       print * ,'instead of the format of the input file: modify your directives'
220       STOP
221 ELSE
222     ! input phase avec obs2mesonh post Oct2005
223   IF (LEN_TRIM(YFILEOUTNAME) == 0) THEN
224     YFILEOUTSUFFIX='obs'
225     YFILEOUTNAME=YFILEGRID
226   ELSE
227     YFILEOUTSUFFIX='NEN'
228     PRINT*,'Output file name=', TRIM(YFILEOUTNAME)
229   ENDIF
230 ENDIF
231 !
232 ! fichier a creer dans WRITEvar
233 YFLAGWRITE='NEW'
234 ! ecriture de ZSBIS dans FILEOUT
235 CALL WRITEVAR(NREADIL,NREADIH,NREADJL,NREADJH,NREADKL,NREADKL,&
236               1,SIZE(XVAR,4),1,SIZE(XVAR,5),1,SIZE(XVAR,6), &
237               'ZSBIS',YFILEOUTNAME,YFLAGWRITE,YFILEOUTSUFFIX,ilocverbia,iret)
238 !  indiquera  a WRITEVAR que le fichier courant en ecriture est deja ouvert
239 YFLAGWRITE='OLD'
240 !-------------------------------------------------------------------------------
241 !
242 !*       3.    LOOP on obs files
243 !              -----------------
244 !    
245 PRINT*,' Loop on the observation files:'
246 DO JOBS=1,10000
247   777 CONTINUE  ! Point de reprise pour le traitement de la 2e composante vent
248   !
249   PRINT*,'- Format of the input observation file:' 
250   PRINT*,'  LL=   n lines Lon,Lat,val '
251   PRINT*,'  ll=   n lines lat,lon,val '
252   PRINT*,'  DLL=  date (YYYYMMDDHHMISS) then n lines Lon,Lat,val '
253   PRINT*,'  Dll=  date (YYYYMMDDHHMISS) then n lines lat,lon,val '
254   PRINT*,'  LLa=  n lines Lon,Lat,alt(m),val'
255   PRINT*,'  lla=  n lines lat,lon,alt(m),val'
256   PRINT*,'  DLLa= date (YYYYMMDDHHMISS) then n lines Lon,Lat,alt(m),val'
257   PRINT*,'  Dlla= date (YYYYMMDDHHMISS) then n lines lat,lon,alt(m),val'
258   PRINT*,'(END to stop)?'
259   READ(5,'(A)',END=88) YLL
260   YLL=ADJUSTL(YLL)
261   CALL WRITEDIR(ILUDIR,YLL)
262   IF (YLL(1:3)=='end' .OR. YLL(1:3)=='END') GO TO 88
263   IF (YLL(1:3)=='LLa') THEN
264     print*, 'format Lon,Lat,alt(m),value'
265     galtobs=.true.
266   ELSE IF (YLL(1:3)=='lla') THEN
267     print*, 'format lat,lon,alt(m),value'
268     galtobs=.true.
269   ELSE IF (YLL(1:4)=='DLLa') THEN
270     print*, 'format Date then Lon,Lat,alt(m),valeur'
271     galtobs=.true.
272   ELSE IF (YLL(1:4)=='Dlla') THEN
273     print*, 'format Date then lat,lon,alt(m),valeur'
274     galtobs=.true.
275   ELSE IF (YLL(1:2)=='LL') THEN
276     print*, 'format Lon,Lat,valeur'
277     galtobs=.false.
278   ELSE IF (YLL(1:2)=='ll') THEN
279     print*, 'format lat,lon,valeur'
280     galtobs=.false.
281   ELSE IF (YLL(1:3)=='DLL') THEN
282     print*, 'format Date then Lon,Lat,value'
283     galtobs=.false.
284   ELSE IF (YLL(1:3)=='Dll') THEN
285     print*, 'format Date then lat,lon,value'    
286     galtobs=.false.
287   ELSE
288     print*, '** incorrect format ',YLL
289     CYCLE
290   ENDIF
291   PRINT*,'- Name of the input observation file ?'
292   READ(5,'(A)',END=88) YFILEOBS
293   YFILEOBS=ADJUSTL(YFILEOBS)
294   CALL WRITEDIR(ILUDIR,YFILEOBS)
295   !
296   !*       3.1   Lecture du fichier d obs a traiter
297   !              ----------------------
298   PRINT*, '- Name of the new field to be created:'
299   PRINT*, '(if it is a wind field you have to name the field : '
300   PRINT*, ' WTxx: the field is localised at vertical flux points, ',&
301           'otherwise at mass points (example : WT10) '
302   PRINT*, ' UTxx: the field (U-component for zonal) will be converted to ',&
303           'MesoNH wind components'
304   PRINT*, 'the V-component must be provided immediately after with VTxx'
305   PRINT*, '?'
306   READ(5,'(A9)',END=88) CGROUP
307   CGROUP=ADJUSTL(CGROUP)
308   CALL WRITEDIR(ILUDIR,CGROUP)
309   CALL LOW2UP(CGROUP)      
310   PRINT*, '- Unit of the new field ?'
311   READ(5,'(A)') YUNITE
312   YUNITE=ADJUSTL(YUNITE)
313   CALL WRITEDIR(ILUDIR,YUNITE)
314   PRINT*, '- Profil of the new field :'
315   PRINT*, '   3D=XYZ '
316   PRINT*, '   2D=XY  (obs altitudes not taken into account)'
317   PRINT*, '   1D=Z   (vertical profil (_PV_ for diaprog) localised at ',&
318            'lat-lon of the 1st obs'
319   PRINT*, ' 1D/2D/3D ?'        
320   READ(5,'(A)') YSTOCK      
321   YSTOCK=ADJUSTL(YSTOCK)
322   CALL WRITEDIR(ILUDIR,YSTOCK)
323   IF ( (YSTOCK == '3D' .OR. YSTOCK == '1D') .AND. .NOT.(galtobs) ) THEN
324       print * ,'** It is not possible to store ',TRIM(YSTOCK),' profil ',&
325                'because no altitude was provided in the input obs file'
326       print *, ' change your inputs:', TRIM(YLL)
327       STOP
328   ENDIF
329   ! tableau de stockage des valeurs des obs et compteur de ces valeurs stockees
330   IF(ALLOCATED(ZOBSinMNH)) DEALLOCATE(ZOBSinMNH)
331   IF(ALLOCATED(ZALTOBS)) DEALLOCATE(ZALTOBS)
332   IF(ALLOCATED(ICPTOBSinMNH)) DEALLOCATE(ICPTOBSinMNH)
333   ! XVAR = futur tableau a ecrire via writevar
334   IF(ALLOCATED(XVAR)) DEALLOCATE(XVAR)
335   IF ( YSTOCK == '3D' .OR. YSTOCK == '1D' ) THEN
336     ALLOCATE(XVAR( SIZE(XZZ,1),SIZE(XZZ,2),SIZE(XZZ,3),1,1,1))
337     ALLOCATE(ZOBSinMNH( SIZE(XZZ,1),SIZE(XZZ,2),SIZE(XZZ,3)))
338     ALLOCATE(ZALTOBS( SIZE(XZZ,1),SIZE(XZZ,2),SIZE(XZZ,3)))
339     ALLOCATE(ICPTOBSinMNH( SIZE(XZZ,1),SIZE(XZZ,2),SIZE(XZZ,3)))
340   ELSE
341     ALLOCATE(XVAR( SIZE(XZZ,1),SIZE(XZZ,2),1,1,1,1))
342     ALLOCATE(ZOBSinMNH( SIZE(XZZ,1),SIZE(XZZ,2),1))
343     ALLOCATE(ZALTOBS( SIZE(XZZ,1),SIZE(XZZ,2),1))
344     ALLOCATE(ICPTOBSinMNH( SIZE(XZZ,1),SIZE(XZZ,2),1))
345   END IF
346   !
347   PRINT*, 'Mesonh field to be created: ', TRIM(CGROUP),' ',TRIM(YUNITE),' ',TRIM(YSTOCK)
348   ! init de la grille verticale Mesonh suivant le nom de variable
349   SELECT CASE (CGROUP(1:1))
350     CASE ('W') 
351       IGRID=4        ! champ d obs sur la grille W
352       ! init du tableau des altitudes XZZ pour la grille masse
353       CALL COMPCOORD_FORDIACHRO(1)            
354       if (ilocverbia > 0 ) then
355         print *,' after COMPCOORD_FORDIACHRO (mass grid for W field)'
356       endif
357       !                             -----------XZZ(k)   grille masse
358       !                                 x W(k)
359       !                             -----------XZZ(k-1) grille masse
360     CASE default 
361       IGRID=1        ! champ d obs sur la grille de masse 
362       ! init du tableau des altitudes XZZ pour la grille W 
363       CALL COMPCOORD_FORDIACHRO(4)            
364       if (ilocverbia > 0 ) then
365         print *,' after COMPCOORD_FORDIACHRO (W grid for mass field)'
366       endif
367       !                             -----------XZZ(k+1) grille W
368       !                                 x T(k)
369       !                             -----------XZZ(k)   grille W
370   END SELECT
371   !
372   !* 3.2.  pour chaque obs lue, recherche de la maille mesonh
373   !        contenant cette obs : cumul dans cette maille
374   !                              mis a jour du compteur d obs par maille
375   !              ----------------------
376   ! 
377   print *,'YFILEOBS=', TRIM(YFILEOBS)
378   CALL CREATLINK('DIROBS',YFILEOBS,'CREAT',ilocverbia)
379   OPEN (UNIT=8,FILE=TRIM(ADJUSTL(YFILEOBS)),STATUS='OLD',FORM='FORMATTED')
380   !
381   ZOBSINMNH(:,:,:)=0.
382   ICPTOBSinMNH(:,:,:)=0
383   IKMAX=1
384   !
385   inbval=0
386   inbvalrej=0
387   if (ilocverbia >= 2 ) print *, 'before reading of input obs file'
388   IF (YLL(1:1)=='D' ) THEN
389     ! lecture de la date
390     READ (8,'(A14)',ERR=886) YDATEOBS
391     YDATEOBS=ADJUSTL(YDATEOBS)
392     ! verification YDATEOBS est une date
393     IF (YDATEOBS(1:4)<='1900' .OR. YDATEOBS(1:4)>='2020' &
394       .OR. YDATEOBS(5:6)<'01' .OR. YDATEOBS(5:6)>'12' &
395       .OR. YDATEOBS(7:8)<'01' .OR. YDATEOBS(7:8)>'31' &
396       .OR. YDATEOBS(9:10)<'00' .OR. YDATEOBS(9:10)>'23' &
397       .OR. YDATEOBS(11:12)<'00' .OR. YDATEOBS(11:12)>'59' &
398       .OR. YDATEOBS(13:14)<'00' .OR. YDATEOBS(13:14)>'59' ) GO TO 887
399     ! Pourquoi cette init ? voir GJ iret=1
400     ! reinitialisation de XDATIME
401     READ(YDATEOBS,'(i4,5I2)') IAN, IMOIS, IJOUR, IHEUR, IMINU, ISEC
402     XSEC=IHEUR*3600.+IMINU*60.+ISEC
403     DO JILOOP=1,13,4
404       XDATIME(JILOOP,1)=IAN
405       XDATIME(JILOOP+1,1)=IMOIS
406       XDATIME(JILOOP+2,1)=IJOUR
407       XDATIME(JILOOP+3,1)=XSEC
408     END DO
409     if (ilocverbia >= 2 ) print *,'XDATIME initialised to:',XDATIME(1:4,1)
410   ENDIF
411   ! lecture de la position et valeur des observations
412   ! boucle infinie : arret sur fin de fichier
413   ! init de la position à -999
414   IIMNH=-999
415   IJMNH=-999
416 ! modif GJ deduction format  JILOOP=0
417   DO
418 ! modif GJ deduction format    IF ( JILOOP == 0 .AND. YSTOCK =='2D'  ) THEN
419 ! modif GJ deduction format       ! le format de donnees peut ne pas contenir alt 
420 ! modif GJ deduction format       ! lecture du premier enregistrement en caracteres puis decodage
421 ! modif GJ deduction format      if (ilocverbia >= 2 ) print *,'Recherche du nombre de variables dans un enregistrement'
422 ! modif GJ deduction format      READ (8,'(A100)') YLIGNE
423 ! modif GJ deduction format      YLIGNE=ADJUSTL(YLIGNE)
424 ! modif GJ deduction format      il=index(YLIGNE,' ')-1
425 ! modif GJ deduction format      READ(Yligne(1:il),*) ZVAR1
426 ! modif GJ deduction format      YLIGNE=ADJUSTL(YLIGNE(il+1:))
427 ! modif GJ deduction format      il=index(YLIGNE,' ')-1
428 ! modif GJ deduction format      READ(Yligne(1:il),*) ZVAR2
429 ! modif GJ deduction format      IF (TRIM(YLL)=='LL' .or. TRIM(YLL)=='DLL') THEN
430 ! modif GJ deduction format        ZOBSLONlu=ZVAR1
431 ! modif GJ deduction format        ZOBSLATlu=ZVAR2
432 ! modif GJ deduction format      ELSE IF (YLL(1:2)=='ll' .or. YLL(1:3)=='Dll') THEN
433 ! modif GJ deduction format        ZOBSLATlu=ZVAR1
434 ! modif GJ deduction format        ZOBSLONlu=ZVAR2
435 ! modif GJ deduction format      ENDIF
436 ! modif GJ deduction format      YLIGNE=ADJUSTL(YLIGNE(il+1:))
437 ! modif GJ deduction format      il=index(YLIGNE,' ')-1
438 ! modif GJ deduction format      READ(Yligne(1:il),*) ZVAR3
439 ! modif GJ deduction format      IF ( LEN_TRIM(YLIGNE(il+1:)) == 0 ) THEN
440 ! modif GJ deduction format        print *,' Champ 2D avec enregistrement sans altitude'
441 ! modif GJ deduction format        ZOBSALTlu=-999
442 ! modif GJ deduction format        ZVALOBS=ZVAR3
443 ! modif GJ deduction format        YLL(LEN_TRIM(YLL):LEN_TRIM(YLL))='A'
444 ! modif GJ deduction format      ELSE
445 ! modif GJ deduction format        ZOBSALTlu=ZVAR3
446 ! modif GJ deduction format        YLIGNE=ADJUSTL(YLIGNE(il+1:))
447 ! modif GJ deduction format        READ(Yligne,*) ZVALOBS
448 ! modif GJ deduction format      ENDIF
449 ! modif GJ deduction format      if (ilocverbia >= 2 ) THEN
450 ! modif GJ deduction format        print *,'1er enregistrement: lat=',ZOBSLATlu,' lon=',ZOBSLONlu
451 ! modif GJ deduction format        print *,'                    alt=',ZOBSALTlu,' var=',ZVALOBS
452 ! modif GJ deduction format      endif
453 ! modif GJ deduction format    ELSE
454 ! modif GJ deduction format      ! lecture d'un enregistrement
455 ! modif GJ deduction format      IF (TRIM(YLL)=='LL' .or. TRIM(YLL)=='DLL') THEN
456 ! modif GJ deduction format        READ (8,*,ERR=887,END=888) ZOBSLONlu,ZOBSLATlu,ZOBSALTlu,ZVALOBS
457 ! modif GJ deduction format      ELSE IF (YLL(1:2)=='ll' .or. YLL(1:3)=='Dll') THEN
458 ! modif GJ deduction format        READ (8,*,ERR=887,END=888) ZOBSLATlu,ZOBSLONlu,ZOBSALTlu,ZVALOBS
459 ! modif GJ deduction format      ELSE IF (TRIM(YLL)=='LA' .or. TRIM(YLL)=='DLA') THEN
460 ! modif GJ deduction format        READ (8,*,ERR=887,END=888) ZOBSLONlu,ZOBSLATlu,ZVALOBS
461 ! modif GJ deduction format        ZOBSALTlu=-999
462 ! modif GJ deduction format      ELSE IF (YLL(1:2)=='lA' .or. YLL(1:3)=='DlA') THEN
463 ! modif GJ deduction format        READ (8,*,ERR=887,END=888) ZOBSLATlu,ZOBSLONlu,ZVALOBS
464 ! modif GJ deduction format        ZOBSALTlu=-999
465 ! modif GJ deduction format      ELSE
466 ! modif GJ deduction format        print * ,' Format des obs =',YLL(1:4),' valeur incorrecte'
467 ! modif GJ deduction format        print *, 'valeurs possibles: ll ou LL ou Dll ou DLL ou llh ou LLh ou Dllh ou DLLh'
468 ! modif GJ deduction format        STOP
469 ! modif GJ deduction format      ENDIF
470 ! modif GJ deduction format    ENDIF
471 ! modif GJ deduction format    JILOOP=1
472     IF (YLL(1:3)=='LLa' .OR. YLL(1:4)=='DLLa') THEN
473       READ (8,*,END=888) ZOBSLONlu,ZOBSLATlu,ZOBSALTlu,ZVALOBS
474     ELSE IF (YLL(1:3)=='lla'.OR. YLL(1:4)=='Dlla' ) THEN
475       READ (8,*,END=888) ZOBSLATlu,ZOBSLONlu,ZOBSALTlu,ZVALOBS
476     ELSE IF (YLL(1:2)=='LL'.OR.YLL(1:3)=='DLL' ) THEN
477       READ (8,*,END=888) ZOBSLONlu,ZOBSLATlu,ZVALOBS
478       ZOBSALTlu= XSPVAL
479     ELSE IF (YLL(1:2)=='ll' .OR. YLL(1:3)=='Dll' ) THEN
480       READ (8,*,END=888) ZOBSLATlu,ZOBSLONlu,ZVALOBS
481       ZOBSALTlu= XSPVAL
482     ELSE
483       print * ,'** Obs format=',YLL(1:4),' is an incorrect value'
484       print *, 'correct values are: ll or LL or Dll or DLL or lla or LLa or Dlla or DLLa'
485       STOP
486     ENDIF  
487
488     ! recupere les coordonnées de l obs sur le plan conforme
489     IF (YSTOCK == '3D' .OR. YSTOCK == '2D' .OR. &
490                          (YSTOCK == '1D' .AND. IIMNH == -999) ) THEN
491       ! recupere pour chaque obs si 2D ou 3D , pour la premiere obs si 1D
492       !(les 2 premiers arg. doivent etre XXHAT et XYHAT (pas XXX et XXY))
493       !! peu importe en masdev4_6 car plus utilises..
494       !CALL SM_XYHAT(XXHAT,XYHAT,XLATORI,XLONORI, &
495       !! XXHAT,XYHAT supprimes en masdev4_7
496       CALL SM_XYHAT(XLATORI,XLONORI, &
497                      ZOBSLATlu,ZOBSLONlu,ZXOBS,ZYOBS)
498       ! quelle est la maille horizontale mesonh qui contient cette obs ?
499       ! XXHAT(I),XXHAT(I+1) = limites X de la maille I
500       ! XYHAT(J),XYHAT(J+1) = limites Y de la maille J
501       IF ( ZXOBS >= XXHAT(2) .AND. ZXOBS <= XXHAT(NIMAX+2-1) .AND.&
502            ZYOBS >= XYHAT(2) .AND. ZYOBS <= XYHAT(NJMAX+2-1) ) THEN
503         IIMNH=MAX(MIN(COUNT(XXHAT(:)<ZXOBS),NIMAX+2-1),2)
504         IJMNH=MAX(MIN(COUNT(XYHAT(:)<ZYOBS),NJMAX+2-1),2)
505       ELSE
506         print * ,'*** The observation at lat,lon ',ZOBSLATlu,ZOBSLONlu,&
507                  'is out of the Mesonh domain, not treated ***'
508         inbvalrej=inbvalrej+1
509         CYCLE
510       ENDIF
511     ELSE
512      if (ilocverbia >= 2 ) then
513        print *, ' Profil ', YSTOCK, ': following obs are at the same localisation',&
514                 ' as the first one ', ZOBSLATlu,ZOBSLONlu
515        print *, 'i,j=', IIMNH,IJMNH
516      endif
517     ENDIF
518
519     if (ilocverbia >= 3 ) then
520        print *, ZXOBS,IIMNH &
521                ,XXX(IIMNH,IGRID),XXX(IIMNH-1,IGRID),XXHAT(IIMNH),&
522                 XXHAT(IIMNH-1),XXHAT(IIMNH+1)
523        print *, ZYOBS,IJMNH &
524                ,XXY(IJMNH,IGRID), XXY(IJMNH-1,IGRID),XYHAT(IJMNH),&
525                 XYHAT(IJMNH-1),XYHAT(IJMNH+1)
526     endif
527     IF ( YSTOCK == '3D' .OR. YSTOCK == '1D' ) THEN
528       ! quelle est la maille verticale mesonh qui contient cette obs ?
529       ! cas des obs a localiser sur la grille de masse
530       ! XZZ_W (K) , XZZ_W(K+1) = limites Z de la maille_masse K
531       ! cas des obs à localiser sur la grille de W
532       ! XZZ_masse (K-1) , XZZ_masse(K) = limites Z de la maille_W K
533       IKMNH=MIN(COUNT(XZZ(IIMNH,IJMNH,:)< ZOBSALTlu),NKMAX+2-1)
534       IF ( IKMNH == 0 .AND. ZVALOBS /= XSPVAL ) THEN
535         print *,'obs under the first model level, stored at',&
536              ' k=1', ZOBSLONlu,ZOBSLATlu,ZOBSALTlu,ZVALOBS
537         IKMNH=1
538       ENDIF
539       IF ( IGRID == 4 ) THEN  ! champ d obs sur la grille W
540         IKMNH=IKMNH+1
541       ENDIF
542       ! stocke le niveau max pour minimiser la taille du tableau a ecrire
543       IKMAX=MAX(IKMAX,IKMNH)
544     ELSE
545       IKMNH=1
546     ENDIF
547     ! stockage
548     if (ilocverbia >= 2 ) then
549       IKM1=MAX(IKMNH-1,1)
550       print *, ZOBSLONlu,ZOBSLATlu,ZOBSALTlu,ZVALOBS,IIMNH,IJMNH,IKMNH &
551              , XZZ(IIMNH,IJMNH,IKMNH),XZZ(IIMNH,IJMNH,IKM1)
552     endif
553     IF (ZVALOBS /= XSPVAL ) THEN
554       if (ilocverbia >= 2 ) then
555         print *,'before TO_COMPUTING_UNITS', ZVALOBS
556       endif
557       CALL TO_COMPUTING_UNITS(CGROUP,YUNITE,ZVALOBS) 
558       if (ilocverbia >= 2 ) then
559         print *,'after TO_COMPUTING_UNITS', ZVALOBS
560       endif
561       ! Voir une amelioration en moyenne ponderee avec la distance :
562       ! serait utile pour des mailles tres grandes
563       if (ilocverbia >=3 ) then
564         print *, 'Storage indexes i,j,k=',IIMNH,IJMNH,IKMNH
565       endif
566       ZALTOBS(IIMNH,IJMNH,IKMNH)=ZOBSALTlu
567       ZOBSinMNH(IIMNH,IJMNH,IKMNH)=ZOBSinMNH(IIMNH,IJMNH,IKMNH)+ZVALOBS
568       ICPTOBSinMNH(IIMNH,IJMNH,IKMNH)=ICPTOBSinMNH(IIMNH,IJMNH,IKMNH)+1
569     ENDIF
570     !
571     inbval=inbval+1
572   ENDDO   ! fin de boucle de lecture du fichier d obs
573 GO TO 888
574 886 CONTINUE
575   print *,' *** WARNING: in reading the date in the obs file ***'
576   print *,' not enough rows (4)'
577   GO TO 888
578 887 CONTINUE
579   print *,' *** WARNING: in reading the obs file ***'
580   print *,'             every record must contains 4 values'
581   print *,'             or 3 in 2D'
582 888   CONTINUE
583   !
584   print *, 'End of reading the input obs file'
585   CLOSE (UNIT=8)
586   ! suppression du lien
587   CALL CREATLINK('DIROBS',YFILEOBS,'CLEAN',ilocverbia)
588   print *, 'number of obs taken into account in the model grid= ', inbval
589   print *, 'number of obs out of domain not taken into account= ', inbvalrej
590   !
591   ! mise a indef des mailles MNH non concernées par les obs
592   !
593   WHERE ( ICPTOBSinMNH(:,:,:) == 0)  
594      ZOBSinMNH(:,:,:)= XSPVAL
595      ZALTOBS(:,:,:)= XSPVAL
596   END WHERE
597   print *, 'number of meshes set to indef= ', COUNT(ICPTOBSinMNH(:,:,:) ==0) 
598   print *, 'number of meshes initialised= ', COUNT(ICPTOBSinMNH(:,:,:) > 0) 
599
600   IF ( (COUNT (ICPTOBSinMNH(:,:,:) > 0) ) == 0 ) THEN
601      print *, '**** no observation is localised into the model grid'
602      print *, ' the field is not written in the output diachronic file'
603      CYCLE
604   ENDIF
605   !
606   ! calcul eventuel de la moyenne des obs incluses dans les mailles mesonh
607   WHERE ( ICPTOBSinMNH(:,:,:) > 0) &
608     ZOBSinMNH(:,:,:)=ZOBSinMNH(:,:,:)/ICPTOBSinMNH(:,:,:)
609   print *, 'end of computation of the average on ',&
610     COUNT(ICPTOBSinMNH(:,:,:) >0) , ' meshes'
611   !
612   ! traitement particulier des composantes du vent
613   SELECT CASE (CGROUP(1:1))
614     CASE ('U','V')
615       IF ( .NOT. ALLOCATED (ZVENTSAVE) ) THEN
616         ALLOCATE(ZVENTSAVE( SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3), &
617                             SIZE(XVAR,4),SIZE(XVAR,5),SIZE(XVAR,6)   ))
618         ALLOCATE(ZWORK3D ( SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3) )) 
619         ALLOCATE(ZWORK3D2( SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3) )) 
620         print *, 'Treatment for the wind: storage of the zonal component 1'
621         ZWORK3D(:,:,:)=ZOBSinMNH(:,:,:)
622         print *, '  and treatement of the Obs file for the 2d component'
623         GO TO 777
624       ELSE
625         print *, 'Treatment for the wind: storage of the meridional component 2'
626         ZWORK3D2(:,:,:)=ZOBSinMNH(:,:,:)
627         CALL UV_TO_ZONAL_AND_MERID(ZWORK3D,ZWORK3D2,0,      &
628                                 PZC=ZVENTSAVE(:,:,:,1,1,1), &
629                                 PMC=XVAR(:,:,:,1,1,1)       )
630         
631         print *,' after UV_TO_ZONAL_AND_MERID'
632         DEALLOCATE( ZWORK3D,ZWORK3D2)
633       ENDIF
634       ! Fin traitement particulier des composantes du vent
635     CASE DEFAULT
636       ! init du champ  passe par module a writevar
637        XVAR(:,:,:,1,1,1)=ZOBSinMNH(:,:,:)
638   ENDSELECT
639   ! 
640   ! init des variables passees par module a writevar
641   NGRIDIA(1)=IGRID
642   CTITRE(1)=CGROUP
643   CCOMMENT(1)='from '//ADJUSTL(YFILEOBS)
644   CUNITE(1)=YUNITE
645   !
646   !*      3.3 Ecriture  du tableau XVAR (module MODD_ALLOC_FORDIACHRO)
647   !           --------------------------------------------------
648   zmini=MINVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
649   zmaxi=MAXVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
650   print * ,' After treatment, min,max of the field ',TRIM(CGROUP),'=', zmini,zmaxi
651   print *,' Writing in diachronic format'
652   if (ilocverbia >= 1 ) then
653     print *,'dimensions of XVAR ', SIZE(XVAR,1) , SIZE(XVAR,2), SIZE(XVAR,3)
654   endif
655   !
656   ivarideb=1
657   ivarifin=SIZE(XVAR,1)
658   ivarjdeb=1
659   ivarjfin=SIZE(XVAR,2)
660   ivarkdeb=1
661   if ( IKMAX <= 2 ) THEN
662     ivarkdeb= IKMAX
663   endif
664   ivarkfin=IKMAX
665   ivartinf=1
666   ivartsup=1
667   ivartrajinf=1
668   ivartrajsup=1
669   ivarprocinf=1
670   ivarprocsup=1
671   IF ( YSTOCK == '1D' ) THEN
672     ! tableaux 1D stockés pour permettre un trace diaprog en profil vertical
673     ivarideb=IIMNH
674     ivarifin=IIMNH
675     ivarjdeb=IJMNH
676     ivarjfin=IJMNH
677     print * ,' Storage of 1D profil, position i,j in the grid=',ivarideb,ivarjdeb
678   ENDIF
679   if (ilocverbia >= 2 ) then
680     print *,'before WRITEVAR',' input arguments ',&
681              ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
682              ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
683              TRIM(CGROUP),' ',TRIM(YFILEOUTNAME),' ',TRIM(YFLAGWRITE),' ',&
684              TRIM(YFILEOUTSUFFIX),&
685              ilocverbia,iret
686   endif
687   CALL WRITEVAR(ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
688                 ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup, &
689                 CGROUP,YFILEOUTNAME,YFLAGWRITE,YFILEOUTSUFFIX,ilocverbia,iret)
690
691   print *, ' WRITEVAR, return value for (',TRIM(CGROUP),')= ',iret
692   IF ( iret /= 0 ) THEN 
693     print *,'** Error when writing in the file: ',TRIM(YFILEOUTNAME)
694     STOP
695   ENDIF              
696   !
697   ! traitement eventuel de la 2e composante du vent
698   IF ( ALLOCATED (ZVENTSAVE) ) THEN
699     XVAR(:,:,:,:,:,:)= ZVENTSAVE(:,:,:,:,:,:)
700     CGROUP='U'//CGROUP(2:)
701     CTITRE(1)='U'//CGROUP(2:)
702     if (ilocverbia >= 2 ) then
703      print *,'before WRITEVAR',' input arguments ',&
704              ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
705              ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
706              TRIM(CGROUP),' ',TRIM(YFILEOUTNAME),' ',TRIM(YFLAGWRITE),' ',&
707              TRIM(YFILEOUTSUFFIX),&
708              ilocverbia,iret
709     endif
710     CALL WRITEVAR(ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
711                 ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup, &
712                 CGROUP,YFILEOUTNAME,YFLAGWRITE,YFILEOUTSUFFIX,ilocverbia,iret)
713
714      print *, ' WRITEVAR, return value for (',TRIM(CGROUP),')= ',iret
715      IF ( iret /= 0 ) THEN 
716        print *,'** Error when writing in the file: ',TRIM(YFILEOUTNAME)
717        STOP
718      ENDIF              
719      DEALLOCATE (ZVENTSAVE)
720   ENDIF
721   !
722   IF (YSTOCK == '2D' ) THEN
723    !IF (COUNT(ZALTOBS(:,:,:) /= XSPVAL) /= 0) THEN
724    IF (galtobs) THEN
725     ! stockage egalement de l altitude des obs comme champ diachronique
726     XVAR(:,:,:,1,1,1)=ZALTOBS(:,:,:)          
727     NGRIDIA(1)=1
728     CTITRE(1)='ALT_'//ADJUSTL(CGROUP)
729     CCOMMENT(1)='from '//ADJUSTL(YFILEOBS)
730     CUNITE(1)='m'          
731     if (ilocverbia >= 2 ) then
732       print *,'before WRITEVAR',' input arguments ',&
733            ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
734            ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
735            TRIM(CGROUP),' ',TRIM(YFILEOUTNAME),' ',TRIM(YFLAGWRITE),' ',&
736            TRIM(YFILEOUTSUFFIX),ilocverbia,iret
737     endif
738     CALL WRITEVAR(ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
739               ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup, &
740               CTITRE(1),YFILEOUTNAME,YFLAGWRITE,YFILEOUTSUFFIX,ilocverbia,iret)
741
742     print *, 'WRITEVAR, return value (',TRIM(CTITRE(1)),')= ',iret
743     IF ( iret /= 0 ) THEN 
744       print *,'** Error when writing in the file: ',TRIM(YFILEOUTNAME)
745       STOP
746     ENDIF              
747    ELSE
748      print * , ' No altitudes in the Obs file: no field ALT_'//ADJUSTL(CGROUP)
749    ENDIF
750   ENDIF
751   !
752   ! + stockage du nombre d obs par point de grille comme champ diachronique
753   XVAR(:,:,:,1,1,1)=ICPTOBSinMNH(:,:,:)          
754   NGRIDIA(1)=1
755   CTITRE(1)='N_'//ADJUSTL(CGROUP)
756   CCOMMENT(1)='from '//ADJUSTL(YFILEOBS)
757   CUNITE(1)='count'          
758   if (ilocverbia >= 2 ) then
759     print *,'before WRITEVAR',' input arguments ',&
760          ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
761          ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
762          TRIM(CTITRE(1)),' ',TRIM(YFILEOUTNAME),' ',TRIM(YFLAGWRITE),' ',&
763          TRIM(YFILEOUTSUFFIX),&
764          ilocverbia,iret
765   endif
766   CALL WRITEVAR(ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
767             ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup, &
768             CTITRE(1),YFILEOUTNAME,YFLAGWRITE,YFILEOUTSUFFIX,ilocverbia,iret)
769
770   print *, 'WRITEVAR return value (',TRIM(CTITRE(1)),')= ',iret              
771   IF ( iret /= 0 ) THEN 
772       print *,'** Error when writing in the file: ',TRIM(YFILEOUTNAME)
773       STOP
774   ENDIF            
775   !
776 ENDDO ! boucle fichier obs a traiter
777 !
778 ! Fin de boucle sur les fichiers d obs
779 ! pour clore le traitement meme si la liste des champs est
780 ! incomplete ( non terminee par END)
781 88  CONTINUE
782 YFILEOBS='END'
783 !
784 !---------------------------------------------------------------------------
785 !
786 !*       4.    Fermeture fichiers
787 !              ------------------
788 !
789 IF ( YFILEOBS(1:3) == 'END' ) THEN
790   PRINT*, 'END -> Close the output file'
791   YFLAGWRITE='CLO'
792   ! dans cet appel seul l argument YFLAGWRITE est pris en compte, tous
793   ! les autres arguments sont ignorés
794   CALL WRITEVAR(ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
795                ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup, &
796                CGROUP,YFILEOUTNAME,YFLAGWRITE,YFILEOUTSUFFIX,ilocverbia,iret)
797   print *, 'WRITEVAR, return value=',iret
798   IF ( iret > 0 ) THEN 
799     print *,'** Error when closing the file: ',TRIM(YFILEOUTNAME)
800   ENDIF            
801 ENDIF
802 !
803 !-------------------------------------------------------------------------------
804 !
805 !*       5.    Fin de programme
806 !              ------------------
807 !
808 99  CONTINUE
809 !
810 !   Suppression de tous les liens eventuellement crees
811 YDUMMYFILE=''
812 CALL CREATLINK(' ',YDUMMYFILE,'CLEAN',ilocverbia)
813 PRINT*, 'The file ',TRIM(YLUDIR),' stores all the input directives'
814 PRINT*, ' you must give a new name to use it again'
815 CLOSE(ILUDIR)
816 !
817 IF (iret==0) THEN
818   print *,'================'
819   IF  (YFILEOUTSUFFIX /= 'NEN' ) THEN
820     PRINT*, 'Output files *',TRIM(YFILEOUTSUFFIX), '.lfi are available'  
821   ELSE
822     PRINT*, 'Output file ', TRIM(YFILEOUTNAME), '.lfi is available'  
823   ENDIF
824   PRINT*, ' Use LCOLAREA=T and LSPOT=T in diaprog to plot the fields'
825 ENDIF
826 !
827 END PROGRAM OBS2MESONH