Philippe 23/02/2016: lfi2cdf: modif for JPHEXT/=1
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / EXTRACTDIA / extractdia.f90
1       PROGRAM  EXTRACTDIA
2 !     ###################
3 !
4 !!****  *EXTRACTDIA* -  lecture d'enregistrements dans fichier diachronique,
5 !                         traitement,
6 !                         ecriture (11 types de format de fichier possibles)
7 !! 
8 !!
9 !!    PURPOSE
10 !!    -------
11
12 !
13 !!**  METHOD
14 !!    ------
15 !!      
16 !     Lecture en entree:
17 !       d'une liste de fichiers diachroniques 
18 !       du format de sortie
19 !       d'une liste de champs a traiter pour chaque fichier diachronique
20 !       d'un zoom selon toutes les directions inclu dans le champ a traiter
21 !         ( seul le zoom selon i,j,k est possible pour le format DIAC)
22 !
23 !     Ecriture en sortie:
24 !       d'un fichier  au format fonction de TYPEOUT c.a.d 
25 !         DIAC= type diachro (un seul fichier contenant toutes
26 !                                       les variables selectionnées)
27 !         LLHV= lon lat alt val (un seul fichier contenant toutes
28 !                                       les variables selectionnées) 
29 !         llhv= lat lon alt val (un seul fichier contenant toutes
30 !                                       les variables selectionnées) 
31 !         ll ou LL zv lon lat  niveau Z val
32 !
33 !         ll ou LL pv lon lat  niveau P val
34 !
35 !----------------------
36 ! AJOUT NOVEMBRE 2009:
37 !----------------------
38 !         IJHV= i j alt val (un seul fichier contenant toutes
39 !                                       les variables selectionnées) 
40 !         jihv= j i alt val (un seul fichier contenant toutes
41 !                                       les variables selectionnées) 
42 !         IJ ou ji zv lon lat  niveau Z val
43 !
44 !         IJ ou ji pv lon lat  niveau P val
45 !----------------------
46 !         FREE= format libre a choisir par l utilisateur (un fichier par variable)
47 !         KCDL ou ZCDL ou PCDL= format CDL (à convertir en netcdf via "tonetcdf")
48 !                               (un seul fichier contenant toutes
49 !                                       les variables selectionnées)
50 !           KCDL si les niveaux verticaux sont les niveaux du modele
51 !           ZCDL si les niveaux verticaux sont des niveaux Z=constante donnes au programme
52 !           PCDL si les niveaux verticaux sont des niveaux P=constante donnes au programme
53 !
54 !  pour les formats *CDL,*Z*,*P*, 2 types de grille horizontale sont possibles:
55 !    'CONF' grille reguliere sur le plan de projection (conforme ou cartesien)
56 !    'LALO' grille reguliere en lat-lon
57 !             dans ce cas les composantes du vent sont transformees
58 !             en composantes zonales et méridiennes.
59 ! sauf pour IJPV, IJZV, jipv, jizv :  CONF obligatoire
60 !!
61 !!    EXTERNAL
62 !!    --------
63 !!          FROM_COMPUTING_UNITS: retour aux unites initiales  avant ecriture
64 !!                               = passage inverse a celui realise par
65 !!                                 TO_COMPUTING_UNITS      
66 !!          appele par writevar,writecdl,writellhv 
67 !!              et par extractdia avant l ecriture au format FREE
68 !!    REFERENCE
69 !!    ---------
70 !!
71 !!    AUTHORS
72 !!    -------
73 !!    I. Mallet , N. Asencio, J. Stein
74 !!
75 !!
76 !!    MODIFICATIONS
77 !!    -------------
78 !!      Original    17/03/2003
79 !       call to dd and ff routines
80 !       call to writeLLHV if LLHV
81 !       clean writevar to delete choice LLHV inside this routine
82 !       add PCDL,LLZV,llzv,LLPV,llpv cases
83 !       allow a zoom 0,0,jdeb,jfin or ideb,ifin,0,0 or 0,0,0,0  05/2005
84 !        add ALT 3Dfield if KCDL, add the LAT and LON 3Dfields if CONF and *CDL
85 !       04/11/2009 (G. Tanguy) : add case IJHV,IJZV, IJPV , JIHV, JIZV, JIPV
86 !       29/03/2011 (G. TANGUY) : add case ZGRB PGRB
87 !       11/07/2014 (G. TANGUY) : correction pour les donnees LES de type SSOl
88 !                                (vlev et field ne correspondaient pas suite à
89 !                                mauvais zoom)
90 !       16/12/2014 (G.DELAUTIER) : ajout cas LLAV llav : altitude au dessus du
91 !       sol
92 !       18/02/2015 (G.DELAUTIER) : ajout cas AGRB : altitude au dessus du
93 !       sol
94 !       Avril 2015 (G.DELAUTIER) : ajout CFIXRESOL pour car GRIB +correction
95 !       pour FF10MAX
96 ! -----------------------------------------------------------------------------
97 !
98 !*       0.    DECLARATIONS
99 !              ------------
100 !
101 ! modules MesoNH
102 USE MODD_CONF, ONLY: NVERB,LCARTESIAN
103 USE MODD_PARAMETERS, ONLY: JPHEXT,JPVEXT,XUNDEF,NUNDEF
104 USE MODD_DIM1, ONLY: NIMAX,NJMAX,NKMAX
105 USE MODD_GRID, ONLY: XLATORI,XLONORI
106 USE MODD_GRID1, ONLY: XZS,XZZ,XLAT,XLON,XXHAT,XYHAT
107 USE MODE_GRIDPROJ  ! subroutines SM_XYHAT et SM_LATLON 
108 USE MODI_UV_TO_ZONAL_AND_MERID
109 USE MODI_HOR_INTERP_4PTS
110 USE MODI_ZINTER
111 USE MODI_PINTER                          
112 ! modules DIACHRO
113 USE MODD_FILES_DIACHRO
114 USE MODN_NCAR,  ONLY: XSPVAL  
115 USE MODD_ALLOC_FORDIACHRO, ONLY: XVAR, &         ! XVAR(i,j,k,t,n,p)
116                                  XTRAJZ, &       ! XTRAJZ(k,t,n)
117                                  XDATIME, &      ! XDATIME(16,t)
118                                  CTITRE, CUNITE,&! CTITRE(p),CUNITE(p)
119                                  NGRIDIA, & ! NGRIDIA(p)
120                                  NGRID
121 USE MODD_COORD, ONLY: XXX,XXY,XXZS, & !  XXX(:,1:7), XXY(:,1:7), XXZS(:,:,1:7)
122                       XXDXHAT,XXDYHAT ! XXDXHAT(:,1:7), XXDYHAT(:,1:7)
123 USE MODD_RESOLVCAR, ONLY: CGROUP, NVERBIA, &
124                           NNDIA, NPROCDIA, NBPROCDIA !pour appel a interp_grids
125 USE MODD_TYPE_AND_LH, ONLY: NIL,NIH,NJL,NJH,NKL,NKH,CTYPE,LICP,LJCP
126 ! modules tools
127 USE MODI_CHANGE_A_GRID 
128 USE MODI_LOW2UP 
129 USE MODI_CREATLINK 
130 USE MODI_DD
131 USE MODI_FF
132 USE MODI_WRITEDIR                                 
133 USE MODI_WRITELLHV
134 USE MODI_WRITEGRIB
135 USE MODI_WRITECDL                                 
136 USE MODI_WRITEVAR                                 
137 USE MODI_FROM_COMPUTING_UNITS
138 USE MODD_READLH
139 USE MODI_INI2LALO
140 USE MODI_INT2LALO
141 !                                 
142 IMPLICIT NONE                       
143 !
144 !*       0.1   Local variables declarations
145 !
146 INTEGER           :: I
147 INTEGER           :: ILUDIR,IRESP
148 INTEGER           :: JLOOP,JI,JJ,JK,J5,J6,J4,JA,JGR,ii
149 ! zoom lu pour les 6 dimensions possibles
150 INTEGER           :: iideb,iifin,ijdeb,ijfin,ikdeb,ikfin
151 REAL              :: zideb,zifin,zjdeb,zjfin
152 INTEGER, dimension(2) :: iloc
153 INTEGER           :: itinf,itsup,itrajinf,itrajsup,iprocinf,iprocsup
154 ! zoom recalcule en fonction des dimensions du champ traite
155 INTEGER           :: ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin
156 INTEGER           :: ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup
157 INTEGER           :: ivarzmin,ivarzmax
158 INTEGER           :: inbvertz,IND_VERT,IND_LL,IND_IJ
159 REAL , allocatable, dimension(:,:,:):: ZWORK3D,ZWORK3D2,zffvent,zdirvent
160 REAL , allocatable, dimension(:,:)  :: zwork2d,zwork2d2
161 REAL , allocatable, dimension(:,:)  :: ZLAT,ZLON
162 REAL , allocatable, dimension(:,:)  :: ZDIFFLON,ZDIFFLAT
163 ! pour traiter les champs budget deja zoomes
164 REAL , allocatable, dimension(:,:,:,:,:,:):: ZVARSAVE                                 
165 ! pour l interpolation verticale a P=cst : pinter
166 REAL , allocatable, dimension(:,:,:) :: ZPABS                          
167 ! pour les interpolations verticales a P ou Z=cst
168 REAL , allocatable, dimension(:,:,:) :: ZVARZCST
169 REAL , allocatable, dimension(:) :: zlistevert
170 INTEGER :: ikdebzint ! premier niveau a traiter
171 !  pour l interpolation sur grille reguliere lat lon
172 REAL , allocatable, dimension(:,:) :: ZNEWLAT,ZNEWLON,ZNEWX,ZNEWY
173 REAL              :: ZDELTALAT,ZDELTALON
174 REAL :: zmini,zmaxi
175 INTEGER           :: inetadd ! compteur de champs supp dans le fichier Netcdf
176 INTEGER           :: IFLAGzcst,IGRID
177 INTEGER           :: IDIM1,IDIM2,I1,I2,IZOOMIDEB,IZOOMIFIN,IZOOMJDEB,IZOOMJFIN
178 INTEGER           :: IAN,IMOIS,IJOUR,IHEURE,IMINUTE,ISECONDE
179
180 INTEGER           :: ilocverbia,iret,iret2,iskip,ISAVENGRIDIA,iarg,INDX,IK
181 CHARACTER(LEN=3)  :: YK
182 !                    flag pour initialiser/ne pas initialiser le zoom d
183 !                      d ecriture : 
184 !                      ne pas initialiser quand ajout par le programme
185 !                      des champs ALT LAT LON qui doivent conserver le
186 !                      zoom de l utilisateur
187 INTEGER           :: ino_init_zoom
188 ! **** la taille des variables caracteres contenant les noms
189 !      de fichiers est obligatoirement de 28 ****
190 CHARACTER(LEN=28) :: YFILEIN,YFILEOUT
191 ! **** la longueur du nom ne doit pas depasser 13 car. si le fichier
192 ! contient des groupes a un seul PROCessus, ou 9 si plusieurs PROCessus ****
193 CHARACTER(LEN=13) :: YGROUP,YGROUP_OLD
194 CHARACTER(LEN=20) :: YGROUP_SAVE
195 CHARACTER(LEN=4)  :: YTYPEOUT
196 CHARACTER(LEN=1)  :: YTYPEOUT3
197 CHARACTER(LEN=3)  :: YSUFFIX_file
198 CHARACTER(LEN=250):: YFMTFREE   ! format ecriture des champs si YTYPEOUT='FREE'
199 CHARACTER(LEN=45) :: YFILEOUTFREE ! nom du fichier de sortie si YTYPEOUT='FREE'
200 CHARACTER(LEN=5)  :: YFLAGREADVAR ,YFLAGWRITE
201 CHARACTER(LEN=4)  :: YOUTGRID  ! grille en sortie:
202                   !CONF pour rester dans le plan conforme,
203                   ! (le logiciel graphique devra réaliser la projection)
204                   !LALO pour passer à lat,lon réguliers
205 CHARACTER(LEN=28) :: YDUMMYFILE
206 CHARACTER(LEN=11) :: YLUDIR      !  Name of the dir file
207 REAL   , DIMENSION(:,:)  ,ALLOCATABLE        :: ZX,ZY             
208 ! GRIB
209 INTEGER :: IND_GRB
210 INTEGER :: ICODCOD ! Parameter grib code
211 INTEGER :: ICODLEV ! grib code for the Type of Level 
212 INTEGER :: ICODOLL ! bottom level if layer 
213 INTEGER :: ICODOLH ! level or top of level if layer
214 CHARACTER(LEN=256):: YINPLINE    ! input agregation line read from Namelist
215 LOGICAL :: LVAR2D
216 INTEGER :: ILEVEL2D ! en option : altitude du champ 2D à coder dans le fichier GRIB
217 LOGICAL :: LLEVEL2D 
218 REAL,DIMENSION(4) :: ZLATLON
219 INTEGER,DIMENSION(4) :: ILATLON
220 INTEGER :: INX,INY
221 REAL,DIMENSION(:,:,:),ALLOCATABLE :: ZALT
222 REAL,DIMENSION(:,:,:),ALLOCATABLE :: zlistevert3D
223 INTEGER :: IZLIST
224 CHARACTER(LEN=1) :: CFIXRESOL
225 REAL :: ZDX_GRB,ZDY_GRB,ZCONTROL
226 !-------------------------------------------------------------------------------
227 !
228 !*       1.     INIT
229 !               ----
230 !
231 !
232 inetadd=0  !compteur de champs supp dans le fichier Netcdf
233 !
234 !Prints : 0=mini 1=debug mode in extractdia, readvar and writevar , writecdl, writellhv
235 !                3=debug mode in routines diachro'
236 ! nverbia= controle des prints dans les routines diachro
237 ilocverbia=0
238
239 ! dans mesonh Xundef est utilise =999.
240 ! dans les routines diachro XSPVAL est utilisé                  
241 XSPVAL=XUNDEF                                    
242 !
243 ! ouverture d un fichier dir ou vont s ecrire les entrees clavier
244 YLUDIR='dirextract'
245 CALL FMATTR(YLUDIR,YLUDIR,ILUDIR,IRESP)
246 OPEN(UNIT=ILUDIR,FILE=YLUDIR,FORM='FORMATTED')
247 !
248 ! Possibilite de definir un zoom d ecriture 
249 !  definition locale du zoom pour extractdia et writevar, writecdl, writellhv
250 iideb=0
251 iifin=0
252 ijdeb=0
253 ijfin=0
254 ikdeb=0
255 ikfin=0
256 itinf=0
257 itsup=0
258 itrajinf=0
259 itrajsup=0
260 iprocinf=0
261 iprocsup=0
262 !
263 !-------------------------------------------------------------------------------
264 !
265 !*       2.     INPUT FILE AND FORMAT
266 !               ---------------------
267 !
268 !*       2.1   name of file and output format
269 !              ------------------------------
270 !
271 PRINT*, '- Name of the diachro file (without .lfi) ?'
272 READ(5,'(A28)') YFILEIN
273 CALL WRITEDIR(ILUDIR,YFILEIN)
274 !
275 PRINT*, '- type of the output file ?'
276 PRINT*, '(DIAC/llhv/llzv/llpv/llav/LLHV/LLZV/LLPV/LLAV/IJHV/IJZV/IJPV/jihv/jizv/jipv/FREE/KCDL/ZCDL/PCDL/ZGRB/PGRB/AGRB)'
277 READ(5,'(A4)')YTYPEOUT
278 CALL WRITEDIR(ILUDIR,YTYPEOUT)
279 PRINT*,'the file ',TRIM(YFILEIN),' will be converted in type ',YTYPEOUT
280 !
281 PRINT*, '- Prints : 0=mini 1=debug mode in extractdia'
282 PRINT*, '                  3=debug mode in routines diachro'
283 PRINT*, '?'
284 READ(5,*)ilocverbia
285 CALL WRITEDIR(ILUDIR,ilocverbia)
286 PRINT*, ' output prints= ',ilocverbia
287 if ( ilocverbia > 2) nverbia=ilocverbia   ! verbosity of diachro routines
288 NVERB=ilocverbia                          ! verbosity of mesonh routines
289 !
290 !*       2.2   other parameters
291 !              ----------------
292 !
293 SELECT CASE (YTYPEOUT)                                   
294   CASE('LLHV','llhv','DIAC','FREE','KCDL','ZCDL','PCDL','llzv','LLZV',&
295           &'llpv','LLPV','IJHV','IJZV','IJPV','jihv','jizv','jipv','ZGRB','PGRB','AGRB',&
296           &'llav','LLAV') ! lecture des choix de l utilisateur
297     IF ( YTYPEOUT == 'FREE' ) THEN
298       PRINT*, '- format of writing for fields ? '
299       PRINT*, '    (fortran syntaxe of FMT in WRITE)'
300       PRINT*,'exemple: (10F9.3) or (8F0.3)'
301       PRINT*, '?'
302       READ(5,'(A)') YFMTFREE
303       CALL WRITEDIR(ILUDIR,YFMTFREE)
304       PRINT*, ' format=', TRIM(YFMTFREE)
305     ENDIF
306     ! lecture du zoom
307     IND_VERT= INDEX(YTYPEOUT(1:4),'Z') + INDEX(YTYPEOUT(1:4),'P') + &
308               INDEX(YTYPEOUT(1:4),'z') + INDEX(YTYPEOUT(1:4),'p') + &
309               INDEX(YTYPEOUT(1:4),'a') + INDEX(YTYPEOUT(1:4),'A')
310     IND_LL= INDEX(YTYPEOUT(1:2),'L') + INDEX(YTYPEOUT(1:2),'l') 
311     IND_IJ= INDEX(YTYPEOUT(1:2),'IJ') + INDEX(YTYPEOUT(1:2),'ji') 
312     IND_GRB=INDEX(YTYPEOUT(1:4),'GRB')
313 print*,YTYPEOUT,IND_IJ
314     IF (IND_LL==0 .AND. IND_GRB==0) THEN
315       IF (IND_VERT/=0) THEN
316         ! cas 'ZCDL','PCDL','jizv','jipv','IJZV','IJPV'
317         PRINT*, '- zoom on the 2 first dimensions: '
318         PRINT*, '              ideb,ifin,jdeb,jfin'
319         PRINT*, '0,0,0,0 for the whole physical domain'
320         PRINT*, '-1,-1,-1,-1 for the whole domain'
321         PRINT*, '?'
322         READ(5,*) iideb,iifin,ijdeb,ijfin
323         CALL WRITEDIR(ILUDIR,iideb)
324         CALL WRITEDIR(ILUDIR,iifin)
325         CALL WRITEDIR(ILUDIR,ijdeb)
326         CALL WRITEDIR(ILUDIR,ijfin)
327       ELSE 
328         ! cas 'DIAC','FREE','KCDL','IJHV','jihv'
329         PRINT*, '- zoom on the 3 first dimensions: '
330         PRINT*, '              ideb,ifin,jdeb,jfin,kdeb,kfin'
331         PRINT*, '0,0,0,0,0,0 for the whole physical domain'
332         PRINT*, '-1,-1,-1,-1,-1,-1 for the whole domain'
333         PRINT*, '?'
334         READ(5,*) iideb,iifin,ijdeb,ijfin,ikdeb,ikfin
335         CALL WRITEDIR(ILUDIR,iideb)
336         CALL WRITEDIR(ILUDIR,iifin)
337         CALL WRITEDIR(ILUDIR,ijdeb)
338         CALL WRITEDIR(ILUDIR,ijfin)
339         CALL WRITEDIR(ILUDIR,ikdeb)
340         CALL WRITEDIR(ILUDIR,ikfin)
341       END IF
342     ELSE
343       ! cas 'llzv','LLZV','llpv','LLPV','llhv','LLHV','llav' 'LLAV'
344       !      'ZGRB','PGRB','AGRB'
345       PRINT*, '- zoom on the 2 first directions: '
346       PRINT*, '              lonmin,lonmax,latmin,latmax'
347       PRINT*, '0.,0.,0.,0. for the whole physical domain'
348       PRINT*, '-1.,-1.,-1.,-1. for the whole domain'
349       PRINT*, '?'
350       READ(5,*) zideb,zifin,zjdeb,zjfin
351       CALL WRITEDIR(ILUDIR,zideb)
352       CALL WRITEDIR(ILUDIR,zifin)
353       CALL WRITEDIR(ILUDIR,zjdeb)
354       CALL WRITEDIR(ILUDIR,zjfin)
355       if(zideb==0. .AND. zifin==0.) then
356         iideb=0 ; iifin=0
357       else if(zideb==-1. .AND. zifin==-1.) then
358         iideb=-1 ; iifin=-1
359       else
360         iideb=-2 ; iifin=-2
361       endif
362       if(zjdeb==0. .AND. zjfin==0.) then
363         ijdeb=0 ; ijfin=0
364       else if(zjdeb==-1. .AND. zjfin==-1.) then
365         ijdeb=-1 ; ijfin=-1
366       else
367         ijdeb=-2 ; ijfin=-2
368       endif
369       IF (IND_GRB/=0) THEN
370         PRINT*,'Do you want to fix resolution in x and y ? (y/n)'
371         PRINT*,'(only available with LALO)'
372         READ(5,*) CFIXRESOL
373         CALL WRITEDIR(ILUDIR,CFIXRESOL)
374         IF (CFIXRESOL=='y') THEN
375           PRINT*,'Enter x resolution (in millidegrees)'
376           READ(5,*) ZDX_GRB
377           PRINT*,'Enter y resolution (in millidegrees)'
378           READ(5,*) ZDY_GRB
379           CALL WRITEDIR(ILUDIR,ZDX_GRB)
380           CALL WRITEDIR(ILUDIR,ZDY_GRB)
381         ENDIF
382
383       ENDIF
384       IF (IND_VERT==0) THEN
385         ! cas 'llhv','LLHV'
386         PRINT*, '- zoom on the 3rd dimension: '
387         PRINT*, '                 kdeb,kfin'
388         PRINT*, '0,0 for the whole physical domain'
389         PRINT*, '-1,-1 for the whole domain'
390         PRINT*, '?'
391         READ(5,*) ikdeb,ikfin
392         CALL WRITEDIR(ILUDIR,ikdeb)
393         CALL WRITEDIR(ILUDIR,ikfin)
394       END IF
395     END IF
396     PRINT*, '- zoom on the 3 last dimensions : '
397     PRINT*, '   itinf,itsup,itrajinf,itrajsup,iprocinf,iprocsup'
398     PRINT*, '0,0,0,0,0,0 for the whole last dimensions'
399     PRINT*, '?'
400     READ(5,*) itinf,itsup,itrajinf,itrajsup,iprocinf,iprocsup
401     CALL WRITEDIR(ILUDIR,itinf)
402     CALL WRITEDIR(ILUDIR,itsup)
403     CALL WRITEDIR(ILUDIR,itrajinf)
404     CALL WRITEDIR(ILUDIR,itrajsup)
405     CALL WRITEDIR(ILUDIR,iprocinf)
406     CALL WRITEDIR(ILUDIR,iprocsup)
407     IF ((iideb==-2) .AND. (ijdeb==-2)) THEN
408       PRINT'(A6,4(E10.4,X),2(I4,X),2(I5,X),4(I4,X))', ' zoom= ',zideb,zifin,zjdeb,zjfin,ikdeb,ikfin&
409                       ,itinf,itsup,itrajinf,itrajsup,iprocinf,iprocsup
410     ELSE
411       PRINT'(A6,6(I4,X),2(I5,X),4(I4,X))', ' zoom= ',iideb,iifin,ijdeb,ijfin,ikdeb,ikfin&
412                       ,itinf,itsup,itrajinf,itrajsup,iprocinf,iprocsup
413     END IF
414     IF (IND_VERT/=0) THEN
415       PRINT*, '- Number of vertical levels for ',YTYPEOUT(IND_VERT:IND_VERT),' interpolation ?'
416       READ(5,*) inbvertz
417       CALL WRITEDIR(ILUDIR,inbvertz)
418       PRINT*, '- Ordered list of these levels (in meters or in hPa): exemple 500 1500 ?'
419       allocate (zlistevert(inbvertz))
420       READ(5,*) zlistevert
421       DO JI=1,inbvertz
422         CALL WRITEDIR(ILUDIR,zlistevert(JI))
423       END DO
424       PRINT*, ' interpolation for the following ',YTYPEOUT(IND_VERT:IND_VERT),' levels='
425       PRINT*, zlistevert 
426     ENDIF
427     YOUTGRID='CONF'
428     IF (YTYPEOUT/='DIAC' .AND. YTYPEOUT/='llhv' .AND. YTYPEOUT/='LLHV' .AND.&
429    & IND_IJ==0) THEN
430       PRINT *,'- Fields in regular LAt/LOn grid'
431       PRINT *,'  or    in regular grid on CONFormal plan (native MesoNH grid) ?'
432       PRINT *,'LALO/CONF ?'
433       READ(5,*) YOUTGRID
434       CALL WRITEDIR(ILUDIR,YOUTGRID)
435       PRINT*, ' Output grid= ', YOUTGRID
436       PRINT*, ''
437       YSUFFIX_file=YTYPEOUT(1:2)//YTYPEOUT(4:4)
438       IF ( YTYPEOUT(2:4) == 'CDL') THEN
439         PRINT*, '!!!!!!!! Warning !!!!!!!!'
440         PRINT*, 'For the CDL type, the dimensions are initialised'
441         PRINT*, ' with those of the first field:'
442         PRINT*, 'the values of the 6 dimensions must be the maximum that'
443         PRINT*, ' will be treated '
444         PRINT*, '!!!!!!!! Warning !!!!!!!!'
445         PRINT*, 'For the CDL type, the coordinates must be the same'
446         PRINT*, ' for all fields'
447         PRINT*, '(stored in the output file with LAT/LON/VLEV groups)'
448         PRINT*, '!!!!!!!!'
449       ENDIF
450     ELSE IF (IND_IJ/=0) THEN ! dans le cas des points de grille on prend les
451                              !  coordonnees conformes
452       YOUTGRID='CONF'
453     ENDIF
454   CASE DEFAULT
455     PRINT*, 'Incorrect value for the output type:',YTYPEOUT
456     PRINT*, 'the following ones are currently available :'
457     PRINT*, 'DIAC,LLHV,llhv,FREE,KCDL,ZCDL,PCDL,llzv,LLZV,llpv,LLPV,llav,LLAV'
458     PRINT*, 'IJHV,IJZV,IJPV,jihv,jizv,jipv,ZGRB,PGRB,AGRB'
459     STOP
460 END SELECT
461
462 !*       2.3   init for input file and output file
463 !              -----------------------------------
464 ! in READVAR, input file must be opened before reading
465 YFLAGREADVAR='OPE'
466 ! in WRITE routine, output file is new
467 YFLAGWRITE='NEW'
468
469 !*       2.4   lecture de la pression pour interpolation
470 !              -----------------------------------------
471 IF (INDEX(YTYPEOUT(1:4),'p')/=0 .OR. INDEX(YTYPEOUT(1:4),'P')/=0 )THEN
472   CALL READVAR('PABST',YFILEIN,YFLAGREADVAR,ilocverbia,iret)
473   IF ( iret /= 0 ) then
474     print *, '- PABST not found, name of the pressure variable ? '
475     read *,YGROUP
476     CALL WRITEDIR(ILUDIR,YGROUP)
477     CALL READVAR(YGROUP,YFILEIN,YFLAGREADVAR,ilocverbia,iret)
478     IF ( iret /= 0 ) then
479       print *,' interpolation at P=cst not possible because PABST and ',TRIM(YGROUP),' are not available'
480       STOP
481     ENDIF
482   ENDIF
483   ! stockage de ZPABS utilise par pinter
484   ALLOCATE ( ZPABS(SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3))) 
485   ZPABS(:,:,:)=XVAR(:,:,:,1,1,1) 
486 ENDIF
487 !
488 !-------------------------------------------------------------------------------
489 !
490 !*       3.    LOOP ON GROUPS IN THE FILE
491 !              --------------------------
492 !
493 DO JGR=1,10000
494   !    
495   !*      3.0  preparation pour la lecture du champ suivant
496   !
497   ino_init_zoom=0
498   IF (IND_GRB==0) THEN
499     PRINT*,'- Name of the group in upper case (13 characters max.)'
500     PRINT*,' (ex: THT or DD or FF or DD10 or FF10 or LAT or LON or VLEV)'
501     PRINT*,'(GROUP for the list of groups, END to stop)?'
502     READ(5,'(A13)',END=88) CGROUP
503     CALL WRITEDIR(ILUDIR,CGROUP)
504     CGROUP=ADJUSTL(CGROUP)
505     CALL LOW2UP(CGROUP)
506   ELSE ! CASE ZGRB or PGRB
507     LLEVEL2D=.FALSE.      
508     LVAR2D=.FALSE. 
509     PRINT*,'- Name of the group in upper case (13 characters max.)'
510     PRINT*,' MesoNH field name, grib parameter indicator'
511     PRINT*,' (ex: UT 131, VT 132, GROUP for the list of groups, END to stop)'
512     PRINT*,' optional : you can add FOR 2D FIELDS ONLY the altitude (in meters)'
513     PRINT*,' of the field after  the grib parameter indicator exple : UT10 131 10'
514     READ(5,'(A)') YINPLINE
515     YINPLINE= TRIM(ADJUSTL(YINPLINE))
516     IF (LEN_TRIM(YINPLINE) == 0) CYCLE ! skip blank line
517     CALL WRITEDIR(ILUDIR,YINPLINE)
518     CALL TAB2SPACE(YINPLINE)
519     ! extract field name
520     INDX= INDEX(YINPLINE,' ')
521     CGROUP= YINPLINE(1:INDX-1)
522     IF (CGROUP=='END') GO TO 88
523     ! 
524     IF (CGROUP /='GROUP') THEN
525       ICODLEV=NUNDEF
526       ICODOLH=NUNDEF
527       ICODOLL=NUNDEF
528       YINPLINE= ADJUSTL(YINPLINE(INDX:))
529       INDX= INDEX(YINPLINE,' ')
530       IF (INDX == 1 ) THEN
531         PRINT*, ' Parameter indicator is missing. ',CGROUP,' not treated.'
532         CYCLE
533       END IF
534       READ(YINPLINE(1:INDX-1),*) ICODCOD
535       IF (NVERB>=5) print*, ' Parameter indicator: ',ICODCOD
536       YINPLINE= ADJUSTL(YINPLINE(INDX:))
537       INDX= INDEX(YINPLINE,' ')
538       IF (INDX /= 1 ) THEN
539         READ(YINPLINE(1:INDX-1),*) ILEVEL2D     
540         PRINT*, 'Level found : ',ILEVEL2D
541         PRINT*, 'it will be only used if the field ',CGROUP,' is 2D'
542         LLEVEL2D=.TRUE.
543       END IF
544
545     ENDIF
546   ENDIF
547   IF (CGROUP=='END') GO TO 88
548   ! point de reprise pour forcer l ecriture des champs VLEV,LAT,LON 
549   ! dans les fichiers netcdf
550 77 CONTINUE
551   YGROUP_SAVE=CGROUP(1:13)
552   YK=''
553   INDX=INDEX(CGROUP,'_K_')
554   IF (INDX/=0) THEN
555     CGROUP=YGROUP_SAVE(1:INDX-1)
556     YK(1:3)=YGROUP_SAVE(INDX+3:INDX+5)
557     READ(YK,'(I3)') IK
558   END IF
559   IF (CGROUP(1:5)/='GROUP') &
560     PRINT*,'you asked for the following record: ',TRIM(CGROUP)
561   !
562   !*      3.1  Lecture et initialisation du tableau XVAR
563   !            passé en module MODD_ALLOC_FORDIACHRO
564   !
565   !
566   !      3.1.1 Cas particulier pour le vent
567   !
568   IF ( CGROUP(1:2) == 'UT' .OR. &
569        CGROUP(1:2) == 'VT' .OR. &
570        CGROUP(1:2) == 'DD' .OR. &
571        CGROUP(1:2) == 'FF' .AND. CGROUP(1:7) /= 'FF10MAX'     )  THEN
572     !
573     IF ( (CGROUP(1:2)=='UT'.OR.CGROUP(1:2)=='VT') .AND. &
574           YOUTGRID(1:4) /= 'LALO'                       ) THEN
575       ! Lecture du champ U ou V sans calcul 
576       ! les composantes du vent restent dans le plan conforme
577       CALL READVAR(CGROUP,YFILEIN,YFLAGREADVAR,ilocverbia,iret)
578     ELSE
579       ! Lecture des 2 composantes du vent  : commence par UM...
580       !(stockees dans les tableaux ZWORK3D et ZWORK3D2)
581       ! max 13 car.
582       YGROUP='UT'//CGROUP(3:13) 
583       CALL READVAR(YGROUP,YFILEIN,YFLAGREADVAR,ilocverbia,iret)
584       IF ( iret /= 0 ) then
585         print *,TRIM(CGROUP),': ',TRIM(YGROUP),' not available'
586         ! echec , on tente UM....
587         YGROUP='UM'//CGROUP(3:13)
588         CALL READVAR(YGROUP,YFILEIN,YFLAGREADVAR,ilocverbia,iret2)
589         IF ( iret2 /= 0 ) then
590           print *,'** no processing for ',TRIM(CGROUP), &
591                   ' because UT and ',TRIM(YGROUP),' are not available'
592           CYCLE
593         ENDIF
594       ENDIF
595       ! allocation du tableau de stockage de la 1e composante du vent
596       ALLOCATE(ZVARSAVE(size(XVAR,1),size(XVAR,2),size(XVAR,3),   &
597                         size(XVAR,4),size(XVAR,5),size(XVAR,6))   )
598       ZVARSAVE=XVAR
599       !
600       ! deuxieme composante VT....
601       YGROUP='VT'//CGROUP(3:13)
602       CALL READVAR(YGROUP,YFILEIN,YFLAGREADVAR,ilocverbia,iret)
603       IF ( iret /= 0 ) then
604         print *,TRIM(CGROUP),': ',TRIM(YGROUP),' not available'
605         ! echec , on tente VM....
606         YGROUP='VM'//CGROUP(3:13)
607    CALL READVAR(YGROUP,YFILEIN,YFLAGREADVAR,ilocverbia,iret2)
608         IF ( iret2 /= 0 ) then
609           print *,'** no processing for ',TRIM(CGROUP), &
610                   ' because VT and ',TRIM(YGROUP),' are not available'
611           CYCLE
612         ENDIF
613         iret=iret2
614       ENDIF
615       !
616       ! Calcul de ff
617       IF (CGROUP(1:2) == 'FF' ) THEN
618         IF (LEN(TRIM(CGROUP)) ==2) THEN
619           YGROUP='VENTFF'
620         ELSE IF (LEN(TRIM(CGROUP)) ==3) THEN
621           YGROUP='VENT'//CGROUP(3:3)//'FF'
622         ELSE IF (LEN(TRIM(CGROUP)) ==4) THEN
623           YGROUP='VENT'//CGROUP(3:4)//'FF'
624         ELSE
625           ! 13 car max
626           YGROUP='VENTFF'//CGROUP(3:9)
627         ENDIF
628         ! allocation du tableau de calcul
629         IF (allocated(ZWORK3D)) DEALLOCATE(ZWORK3D)
630         ALLOCATE(ZWORK3D(size(XVAR,1),size(XVAR,2),size(XVAR,3)))
631         ZWORK3D(:,:,:)=XSPVAL
632         DO J6=1,SIZE(XVAR,6)
633           IGRID=NGRIDIA(J6)
634           DO J5=1,SIZE(XVAR,5)
635           DO J4=1,SIZE(XVAR,4)
636             CALL FF (ZVARSAVE(:,:,:,J4,J5,J6),XVAR(:,:,:,J4,J5,J6),ZWORK3D, &
637                      JPVEXT,JPHEXT,IGRID)
638             XVAR(:,:,:,J4,J5,J6)=ZWORK3D(:,:,:)
639           END DO
640           END DO
641           ! initialisation des variables necessaires a l ecriture
642           CGROUP=YGROUP
643           CTITRE(J6)=YGROUP
644           NGRIDIA(J6)=1
645         END DO
646         DEALLOCATE(ZWORK3D)
647         ! Calcul de dd par rapport au Nord geographique
648       ELSE IF (CGROUP(1:2) == 'DD') THEN
649         IF (CTYPE=='CART' .OR. CTYPE=='MASK' .OR. CTYPE=='SPXY') THEN 
650           IF (LEN(TRIM(CGROUP)) ==2) THEN
651             YGROUP='VENTDD'
652           ELSE IF (LEN(TRIM(CGROUP)) ==3) THEN
653             YGROUP='VENT'//CGROUP(3:3)//'DD'
654           ELSE IF (LEN(TRIM(CGROUP)) ==4) THEN
655             YGROUP='VENT'//CGROUP(3:4)//'DD'
656           ELSE
657           ! 13 car max
658             YGROUP='VENTDD'//CGROUP(3:9) 
659           ENDIF
660           ! allocation du tableau de calcul
661           IF (allocated(ZWORK3D)) DEALLOCATE(ZWORK3D)
662           ALLOCATE(ZWORK3D(size(XVAR,1),size(XVAR,2),size(XVAR,3)))
663           DO J6=1,SIZE(XVAR,6)
664             IGRID=NGRIDIA(J6)
665             DO J5=1,SIZE(XVAR,5)
666             DO J4=1,SIZE(XVAR,4)
667               iskip=1 ! tous les points de grille
668               CALL DD(ZVARSAVE(:,:,:,J4,J5,J6),XVAR(:,:,:,J4,J5,J6),ZWORK3D, &
669                       iskip,IGRID,PLON=XLON(NIL:NIH,NJL:NJH))
670               XVAR(:,:,:,J4,J5,J6)=ZWORK3D(:,:,:)
671             END DO
672             END DO
673             ! initialisation des variables necessaires a l ecriture
674             CGROUP=YGROUP
675             CTITRE(J6)=YGROUP
676             CUNITE(J6)='degrees'
677             NGRIDIA(J6)=1
678           END DO
679           DEALLOCATE(ZWORK3D)
680         ELSE
681           print *,'** processing of ',TRIM(CGROUP),' is not performed for CTYPE= ',CTYPE
682           CYCLE
683         ENDIF
684       ELSE IF (CGROUP(1:2) == 'UT' .OR. CGROUP(1:2) == 'VT') THEN
685         IF (CTYPE=='CART' .OR. CTYPE=='MASK' .OR. CTYPE=='SPXY') THEN 
686         ! Calcul des composantes zonale et meridienne
687         !(YOUTGRID(1:4) == 'LALO') avec la routine UV_TO_ZONAL_AND_MERID
688           print*,' Translate to meridional and zonal wind components'
689           ALLOCATE(ZWORK3D(size(XVAR,1),size(XVAR,2),size(XVAR,3)))
690           ALLOCATE(ZWORK3D2(size(XVAR,1),size(XVAR,2),size(XVAR,3)))
691           IF (ilocverbia >= 3 ) then
692             print *,'before UV_TO_ZONAL_AND_MERID KGRID=23'
693             print'(A31,3(I5,X))',' dimensions of the input arrays',size(ZVARSAVE,1),&
694                                       size(ZVARSAVE,2),size(ZVARSAVE,3)
695             print'(3(I5,X))',size(XVAR,1),size(XVAR,2),size(XVAR,3)
696             print'(A32,3(I5,X))',' dimensions of the output arrays',size(ZWORK3D,1),&
697                                        size(ZWORK3D,2),size(ZWORK3D,3)
698             print'(3(I5,X))',size(ZWORK3D2,1),size(ZWORK3D2,2),size(ZWORK3D2,3)
699           ENDIF
700           DO J6=1,SIZE(XVAR,6)
701             DO J5=1,SIZE(XVAR,5)
702             DO J4=1,SIZE(XVAR,4)
703               CALL UV_TO_ZONAL_AND_MERID(ZVARSAVE(:,:,:,J4,J5,J6), &
704                                          XVAR(:,:,:,J4,J5,J6),     &
705                                          23,PZC=ZWORK3D,PMC=ZWORK3D2)
706               IF (CGROUP(1:1) == 'U' ) THEN
707                 XVAR(:,:,:,J4,J5,J6)=ZWORK3D(:,:,:)
708               ENDIF
709               IF (CGROUP(1:1) == 'V' ) THEN
710                 XVAR(:,:,:,J4,J5,J6)=ZWORK3D2(:,:,:)
711               ENDIF
712             END DO
713             END DO
714           END DO
715           IF (ilocverbia >= 3 ) then
716             print *,'after UV_TO_ZONAL_AND_MERID KGRID=23'
717           END IF
718           ! Stockage dans le tableau XVAR qui est le tableau ecrit
719           ! de la composante souhaitée
720           IF (CGROUP(1:1) == 'U' ) THEN
721             print *, ' U zonal wind component'
722             IF (LEN(TRIM(CGROUP)) ==2) THEN
723               YGROUP='UZON'
724             ELSE IF (LEN(TRIM(CGROUP)) ==3) THEN
725               YGROUP='U'//CGROUP(3:3)//'ZON'
726             ELSE IF (LEN(TRIM(CGROUP)) ==4) THEN
727               YGROUP='U'//CGROUP(3:4)//'ZON'
728             ELSE
729               ! 13 car max
730               YGROUP='UZON'//CGROUP(3:9)
731             ENDIF
732             CTITRE(:)='U zonal wind component'
733           ELSE IF (CGROUP(1:1) == 'V' ) THEN
734             print *, ' V meridian wind component'
735             IF (LEN(TRIM(CGROUP)) ==2) THEN
736               YGROUP='VMED'
737             ELSE IF (LEN(TRIM(CGROUP)) ==3) THEN
738               YGROUP='V'//CGROUP(3:3)//'MED'  
739             ELSE IF (LEN(TRIM(CGROUP)) ==4) THEN
740               YGROUP='V'//CGROUP(3:4)//'MED'
741             ELSE
742               ! 13 car max
743               YGROUP='VZON'//CGROUP(3:9)
744             END IF
745             CTITRE(:)='V meridian wind component'
746           ENDIF
747           CGROUP=YGROUP
748           NGRIDIA(:)=1  ! UZON et VMED en grille de masse
749           DEALLOCATE(ZWORK3D,ZWORK3D2)
750         ELSE
751           print *,' No processing of UZON and VMED for CTYPE= ',CTYPE
752           CYCLE
753         ENDIF
754       ENDIF
755       DEALLOCATE(ZVARSAVE)
756     ENDIF
757   !
758   !      3.1.2 LATitude ou LONgitude de chaque point de la grille conforme
759   !
760   ELSE IF (CGROUP(1:3)=='LAT' .OR. CGROUP(1:3)=='LON') THEN
761     print *, 'LAT/LON asked and YFLAGREADVAR=', YFLAGREADVAR
762    IF ( YFLAGREADVAR /= 'NOP') THEN
763     ! Lecture d un champ 2D quelconque pour initialiser XLAT et XLON
764     CALL READVAR('ZSBIS',YFILEIN,YFLAGREADVAR,ilocverbia,iret)
765     IF ( iret /= 0 ) then
766      ! cas de fichier diachronique sans ZSBIS
767       print *, '- Name of one group in upper case '
768       read *,YGROUP
769       CALL WRITEDIR(ILUDIR,YGROUP)
770       CALL LOW2UP(YGROUP)
771       CALL READVAR(YGROUP,YFILEIN,YFLAGREADVAR,ilocverbia,iret)
772       IF ( iret /= 0 ) then
773          print * ,'**group ', TRIM(YGROUP) , 'not found'
774          stop
775       ENDIF
776     ENDIF
777    ENDIF
778     ! init du tableau XVAR au champ souhaite
779     DEALLOCATE(XVAR)
780     ALLOCATE(XVAR(size(XLAT,1),size(XLAT,2),1,1,1,1) )
781     IF (CGROUP(1:3)=='LAT') THEN
782       XVAR(:,:,1,1,1,1)=XLAT(:,:)
783       CTITRE(1)='latitudes'
784       CUNITE(1)='degrees_north'
785     ELSE IF (CGROUP(1:3)=='LON') THEN
786       XVAR(:,:,1,1,1,1)=XLON(:,:)
787       CTITRE(1)='longitudes'
788       CUNITE(1)='degrees_east'
789     ENDIF
790   !
791   !      3.1.3 Altitude de chaque point de la grille conforme
792   !
793   ELSE IF (CGROUP(1:4)=='VLEV') THEN
794     print *, 'VLEV asked and YFLAGREADVAR=', YFLAGREADVAR
795     IF(CTYPE=='SSOL'.OR.CTYPE=='DRST'.OR.CTYPE=='RAPL'.OR.CTYPE=='RSPL') THEN
796       IF ( YFLAGREADVAR == 'NOP') THEN
797       ! altitude des niveaux du groupe precedent dans XTRAJZ
798         print *,'warning, for CTYPE=',CTYPE,' Vertical LEVels of previous group (',TRIM(YGROUP_OLD),')'
799         DEALLOCATE(XVAR)
800         ALLOCATE(XVAR(1,1,size(XTRAJZ,1),1,1,1))
801         XVAR(1,1,:,1,1,1)=XTRAJZ(:,1,1)
802       ELSE
803         print*,'** no processing with VLEV at the first group'
804         GOTO 99
805       ENDIF
806     ELSE
807       IF ( YFLAGREADVAR /= 'NOP') THEN
808         ! Lecture d un champ 2D quelconque pour initialiser les tableaux XZZ
809         CALL READVAR('ZSBIS',YFILEIN,YFLAGREADVAR,ilocverbia,iret)
810         IF ( iret /= 0 ) then
811           ! cas de fichier diachronique sans ZSBIS
812           print *, '- Name of one group in upper case '
813           read *,YGROUP
814           CALL WRITEDIR(ILUDIR,YGROUP)
815           CALL LOW2UP(YGROUP)
816           CALL READVAR(YGROUP,YFILEIN,YFLAGREADVAR,ilocverbia,iret)
817           IF ( iret /= 0 ) then
818             print * ,'** group ', TRIM(YGROUP) , 'not found'
819             stop
820           ENDIF
821         ENDIF
822       ENDIF
823       ! init de XZZ a la grille de masse ( par defaut readvar 
824       ! l initialise a la grille 4 des  vitesse verticales W)
825       CALL COMPCOORD_FORDIACHRO(1)
826       ! init du tableau XVAR au champ souhaite
827       DEALLOCATE(XVAR)
828       ALLOCATE(XVAR(size(XZZ,1),size(XZZ,2),size(XZZ,3),1,1,1))
829       XVAR(:,:,:,1,1,1)=XZZ(:,:,:)
830       ! retour au XZZ grille 4
831       CALL COMPCOORD_FORDIACHRO(4)
832     ENDIF
833     CTITRE(1)='model levels altitudes ASL'
834     CUNITE(1)='meters'
835   !
836   !      3.1.4 Default case
837   !
838   ELSE
839     !
840     ! Lecture du  champ CGROUP et stockage dans XVAR
841     ! + Initialisation (si YFLAGREADVAR='OPE') des variables
842     ! des modules (cf USE en debut de programme)
843     ! Appel a menu_diachro pour la liste des groupes si CGROUP(1:5)=='GROUP'
844     CALL READVAR(CGROUP,YFILEIN,YFLAGREADVAR,ilocverbia,iret)
845     IF (CGROUP(1:5)=='GROUP') CYCLE
846     !
847   ENDIF
848   !
849   IF ( iret == 0 ) THEN
850     zmini=MINVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
851     zmaxi=MAXVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
852     print * ,' After read, min,max of the variable ',TRIM(CGROUP),'=', zmini,zmaxi
853     !       
854     !*      3.2  Init des bornes min max du zoom en fonction des
855     !            dimensions du tableau XVAR traite
856     !
857    IF ( ino_init_zoom == 0) THEN
858     IF (iideb == 0 .AND. iifin == 0 ) THEN
859       ivarideb=NREADIL ; ivarifin=NREADIH
860       IF (ivarideb/=ivarifin) THEN  ! domI/=1
861         ivarideb=MAX(1+JPHEXT,NREADIL) 
862         ivarifin=MIN(SIZE(XVAR,1)-JPHEXT,NREADIH)
863       ENDIF
864     ELSE IF (iideb == -1 .AND. iifin == -1 ) THEN
865       ivarideb=MAX(1,NREADIL) 
866       ivarifin=MIN(SIZE(XVAR,1),NREADIH)
867     ELSE IF (iideb == -2 .AND. iifin == -2 ) THEN
868       ivarideb=-2
869       iideb=1+JPHEXT
870       IF (zideb >= minval(XLON)) THEN
871         DO JJ=1,SIZE(XLON,2)
872           ivarideb=MAX(MIN(COUNT(XLON(:,JJ)<zideb),SIZE(XLON,1)),iideb)
873           iideb=ivarideb
874         END DO
875       ENDIF
876       ivarifin=-2
877       iifin=1+JPHEXT
878       IF (zifin <= maxval(XLON)) THEN
879         DO JJ=1,SIZE(XLON,2)
880           ivarifin=MAX(MIN(COUNT(XLON(:,JJ)<zifin),SIZE(XLON,1)),iifin)
881           iifin=ivarifin
882         END DO
883       ENDIF
884     ELSE
885       ivarideb=max(iideb,NREADIL)
886       ivarifin=min(iifin,NREADIH)
887       ivarideb=min(ivarideb,ivarifin)
888     ENDIF
889     IF(ijdeb == 0 .AND. ijfin == 0) THEN
890       ivarjdeb=NREADJL ; ivarjfin=NREADJH
891       IF (ivarjdeb/=ivarjfin) THEN  ! domJ/=1
892         ivarjdeb=MAX(1+JPHEXT,NREADJL)
893         ivarjfin=MIN(SIZE(XVAR,2)-JPHEXT,NREADJH)
894       ENDIF
895     ELSE IF (ijdeb == -1 .AND. ijfin == -1 ) THEN
896       ivarjdeb=MAX(1,NREADJL)
897       ivarjfin=MIN(SIZE(XVAR,2),NREADJH)
898     ELSE IF (ijdeb == -2 .AND. ijfin == -2 ) THEN
899       ivarjdeb=-2
900       ijdeb=1+JPHEXT
901       IF (zjdeb >= minval(XLAT)) THEN
902         DO JI=1,SIZE(XLAT,1)
903           ivarjdeb=MAX(MIN(COUNT(XLAT(JI,:)<zjdeb),SIZE(XLAT,2)),ijdeb)
904           ijdeb=ivarjdeb
905         END DO
906       ENDIF
907       ivarjfin=-2
908       ijfin=1+JPHEXT
909       IF (zjfin <= maxval(XLAT)) THEN
910         DO JI=1,SIZE(XLAT,1)
911           ivarjfin=MAX(MIN(COUNT(XLAT(JI,:)<zjfin),SIZE(XLAT,2)),ijfin)
912           ijfin=ivarjfin
913         END DO
914       ENDIF
915     ELSE
916       ivarjdeb=max(ijdeb,NREADJL)
917       ivarjfin=min(ijfin,NREADJH)
918       ivarjdeb=min(ivarjdeb,ivarjfin)
919     ENDIF
920     IF(ivarideb==-2 .OR. ivarifin==-2 .OR. ivarjdeb==-2 .OR.  ivarjfin==-2) THEN
921       print *,'****zoom provided is not included in the FM-file grid'
922       print *,'LON (zoom: ',zideb,zifin,') (file: ',minval(XLON),maxval(XLON)
923       print *,'LAT (zoom: ',zjdeb,zjfin,') (file: ',minval(XLAT),maxval(XLAT)
924       GOTO 99
925     ENDIF
926     IF (IND_VERT/=0) THEN
927       ivarzmin=1   ; ivarzmax=inbvertz
928     ELSE
929       ivarzmin=MAX(1,NREADKL)  ; ivarzmax=MIN(SIZE(XVAR,3),NREADKH)
930       inbvertz=ivarzmax-ivarzmin+1
931     ENDIF
932     IF (ikdeb == 0 .AND. ikfin == 0 ) THEN
933       ivarkdeb=NREADKL ; ivarkfin=NREADKH
934       IF (ivarkdeb/=ivarkfin .AND. CTYPE/='SSOL') THEN  ! domK/=1
935         ivarkdeb=MAX(1+JPVEXT,NREADKL)
936         ivarkfin=min(ivarzmax,SIZE(XVAR,3)-JPVEXT)
937       ENDIF
938     ELSEIF (ikdeb == -1 .AND. ikfin ==-1 ) THEN
939       ivarkdeb=ivarzmin
940       ivarkfin=ivarzmax
941     ELSE
942       ivarkdeb=max(ikdeb,ivarzmin)
943       ivarkfin=min(ikfin,ivarzmax)
944       ivarkdeb=min(ivarkdeb,ivarkfin)
945     ENDIF   
946     IF (INDX/=0) THEN
947       ivarkdeb=IK ; ivarkfin=IK
948     END IF
949    ENDIF
950
951     IF (itinf == 0 .AND. itsup == 0 ) THEN
952       ivartinf=1 ; ivartsup=SIZE(XVAR,4)
953     ELSE
954       ivartinf=max(itinf,1)
955       ivartsup=min(itsup,SIZE(XVAR,4))
956       ivartinf=min(ivartinf,ivartsup)
957     ENDIF
958     IF (itrajinf == 0 .AND. itrajsup == 0 ) THEN
959       ivartrajinf=1 ; ivartrajsup=SIZE(XVAR,5)
960     ELSE
961       ivartrajinf=max(itrajinf,1)
962       ivartrajsup=min(itrajsup,SIZE(XVAR,5))
963       ivartrajinf=min(ivartrajinf,ivartrajsup)
964     ENDIF
965     IF (iprocinf == 0 .AND. iprocsup == 0 ) THEN
966       ivarprocinf=1 ; ivarprocsup=SIZE(XVAR,6)
967     ELSE
968       ivarprocinf=max(iprocinf,1)
969       ivarprocsup=min(iprocsup,SIZE(XVAR,6))
970       ivarprocinf=min(ivarprocinf,ivarprocsup)
971     ENDIF
972     if (ilocverbia > 0 ) then
973       PRINT*,' Zoom limits initialized with:'
974       PRINT'(A53,6(I4,X))','ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin',&
975             ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin 
976       PRINT'(A53,6(I4,X))','ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocfin',&
977             ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup 
978     endif
979     !
980     !*      3.3  Ecriture  du tableau XVAR (module MODD_ALLOC_FORDIACHRO) 
981     !
982     print *,' Write with the format ', YTYPEOUT(1:4)
983     SELECT CASE(YTYPEOUT(1:4))
984       !
985       CASE('DIAC')
986         CALL WRITEVAR(ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
987                       ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,  &
988                       CGROUP,YFILEIN,YFLAGWRITE,'2  ',ilocverbia,iret)
989         if (ilocverbia > 0 ) then
990           print'(A17,I2))','WRITEVAR return= ',iret
991         end if
992       !
993       CASE('FREE')
994         if (ilocverbia >= 0 ) then
995           print*,' format ',YTYPEOUT
996           print'(A53,X,A50,6(I4,X),2(I6,X),4(I4,X))',&
997                  ' domaine for writting : ideb,ifin,jdeb,jfin,kdeb,kfin', &
998                  ',itinf,itsup,itrajinf,itrajsup,iprocinf,iprocsup= ', &
999               ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
1000               ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup 
1001         endif
1002         !  Retour aux unites initiales si necessaire
1003         CALL FROM_COMPUTING_UNITS(CGROUP,CUNITE(1)) 
1004         !
1005         YFILEOUTFREE=ADJUSTL(ADJUSTR(YFILEIN)//'.'//ADJUSTL(ADJUSTR(CGROUP)))
1006         OPEN (UNIT=7,STATUS='NEW',FORM='FORMATTED',FILE=YFILEOUTFREE)
1007         ! a. Ecriture de l entete
1008         !temps courant
1009         IAN=XDATIME(13,1)
1010         IMOIS=XDATIME(14,1)
1011         IJOUR=XDATIME(15,1)
1012         IHEURE=XDATIME(16,1)/3600
1013         IMINUTE=(XDATIME(16,1)-(IHEURE*3600))/60
1014         ISECONDE=ISECONDE-(IHEURE*3600)-(IMINUTE*60)
1015         WRITE(7,FMT='(6(I4,X),2(I6,X),4(I4,X),4(I4,X),A42,A33)') ivarideb,&
1016                    ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
1017                    ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1018                    IAN,IMOIS,IJOUR,IHEURE,IMINUTE ,&
1019                    'format ligne1=  12 Indices (.deb .fin) du ',&
1020                    'tableau  an mois jour hUTC minute'
1021         ! b. ecriture des données au fmt choisi par l utilisateur
1022         WRITE(7,FMT=YFMTFREE) &
1023          XVAR(ivarideb:ivarifin,ivarjdeb:ivarjfin,ivarkdeb:ivarkfin,&
1024               ivartinf:ivartsup,ivartrajinf:ivartrajsup,ivarprocinf:ivarprocsup)
1025         PRINT*,'File ',TRIM(YFILEOUTFREE),' available'
1026         CLOSE(7)
1027       !
1028       CASE('LLHV','llhv','IJHV','jihv')
1029         IF (CTYPE == 'SSOL') THEN
1030           ALLOCATE(ZALT(1,1,SIZE(XTRAJZ,1)))
1031           ZALT(1,1,:)=XTRAJZ(:,1,1)
1032           CALL WRITELLHV(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1033                        ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1034                        ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1035                        CGROUP,YFILEIN,YFLAGWRITE,YTYPEOUT,&
1036                        ilocverbia,iret,PALT=ZALT)
1037         ELSE
1038           CALL WRITELLHV(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1039                        ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1040                        ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1041                        CGROUP,YFILEIN,YFLAGWRITE,YTYPEOUT,&
1042                        ilocverbia,iret)     
1043         ENDIF  
1044         if (ilocverbia > 0 ) then
1045           print*,' WRITELLHV return= ',iret
1046         end if
1047       !
1048       CASE('KCDL','ZCDL','PCDL','LLZV','LLPV','llpv','llzv',&
1049              & 'IJZV','jizv','IJPV','jipv','llav','LLAV')
1050         ! replace field at mass points
1051         IF ( CGROUP /= 'VLEV' ) THEN
1052           If (ALLOCATED(ZWORK3D))DEALLOCATE(ZWORK3D)
1053           If (ALLOCATED(ZWORK3D2))DEALLOCATE(ZWORK3D2)
1054           ALLOCATE(ZWORK3D(SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3)))
1055           ALLOCATE(ZWORK3D2(SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3)))
1056           DO J6=ivarprocinf,ivarprocsup
1057             IGRID=NGRIDIA(J6)
1058             IF(SIZE(XVAR,3)/=1 .OR. IGRID/=4) THEN 
1059               ! pas d interpolation verticale pour champ 2D
1060               DO J5=ivartrajinf,ivartrajsup
1061                 DO J4=ivartinf,ivartsup
1062                   ZWORK3D(:,:,:)=XVAR(:,:,:,J4,J5,J6)
1063                   print'(A29,3(X,I4))',' mass point grid for J4,J5,J6=',J4,J5,J6
1064                   CALL CHANGE_A_GRID(ZWORK3D,IGRID,ZWORK3D2)
1065                   NGRIDIA(J6)=IGRID
1066                   ! IGRID=1 en sortie de change_a_grid
1067                   XVAR(:,:,:,J4,J5,J6)=ZWORK3D2(:,:,:)               
1068                 ENDDO
1069               ENDDO
1070             ENDIF
1071           ENDDO
1072           DEALLOCATE(ZWORK3D,ZWORK3D2)
1073         ENDIF
1074         !
1075         ! a. reinit avant ecriture de la grille verticale correspondant a la
1076         !grille de masse sur laquelle le champ a ete interpole
1077         IFLAGzcst=0
1078         IF (IND_VERT/=0) THEN
1079           IF ( CGROUP == 'VLEV' ) THEN
1080             ! ecriture de la liste des niveaux verticaux 
1081             IFLAGzcst=1
1082             DEALLOCATE(XVAR)
1083             allocate(XVAR(1,1,inbvertz,1,1,1))
1084             XVAR(1,1,:,1,1,1)=zlistevert
1085             ivarideb=1 ; ivarifin=1
1086             ivarjdeb=1 ; ivarjfin=1
1087             ivarkdeb=1 ; ivarkfin=inbvertz
1088             CTITRE(1)='vertical_levels'
1089             CUNITE(1)='user choice'
1090             IF ( YTYPEOUT(IND_VERT:IND_VERT) == 'z' .OR.  YTYPEOUT(IND_VERT:IND_VERT) == 'Z' ) THEN
1091                 CUNITE(1)='km'
1092                 XVAR=XVAR*0.001
1093             ENDIF
1094             IF ( YTYPEOUT(IND_VERT:IND_VERT) == 'p' .OR.  YTYPEOUT(IND_VERT:IND_VERT) == 'P' ) THEN
1095               CUNITE(1)='hPa'
1096             ENDIF
1097           ENDIF
1098         ! b. interpolation eventuelle selon la verticale
1099           IF( SIZE(XVAR,3)>1 .AND. CGROUP /= 'VLEV' ) THEN
1100             ! VLEV, LON, LAT et chps 2D ne passent pas cette partie 
1101             if (ilocverbia >= 0 ) then
1102               print*,' Interpolations on ',inbvertz,' ', &
1103                      YTYPEOUT(IND_VERT:IND_VERT),'-levels'
1104             endif
1105             if (ilocverbia >= 1 .AND. IND_VERT/=0) THEN
1106               print*,'levels= ',zlistevert 
1107             endif
1108             ALLOCATE(ZVARSAVE(size(XVAR,1),size(XVAR,2),size(XVAR,3),   &
1109                               size(XVAR,4),size(XVAR,5),size(XVAR,6))   )
1110             ZVARSAVE=XVAR
1111             ALLOCATE(ZWORK3D(SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3)))
1112             ALLOCATE(ZVARZCST(SIZE(XVAR,1),SIZE(XVAR,2),inbvertz))
1113             DEALLOCATE(XVAR)
1114             ALLOCATE(XVAR(SIZE(ZVARSAVE,1),SIZE(ZVARSAVE,2),SIZE(ZVARZCST,3),&
1115                           size(ZVARSAVE,4),size(ZVARSAVE,5),size(ZVARSAVE,6)))
1116             DO J6=ivarprocinf,ivarprocsup
1117               IGRID=NGRIDIA(J6)
1118               ! init du tableau des altitudes  XZZ pour la grille= IGRID
1119               CALL COMPCOORD_FORDIACHRO(IGRID)
1120               DO J5=ivartrajinf,ivartrajsup
1121                 DO J4=ivartinf,ivartsup
1122                   ZWORK3D(:,:,:)=ZVARSAVE(:,:,:,J4,J5,J6)
1123                   ikdebzint=2
1124                   IF (INDEX(YTYPEOUT(1:4),'Z')/=0 .OR. INDEX(YTYPEOUT(1:4),'z')/=0) THEN
1125                     CALL ZINTER(ZWORK3D,XZZ,ZVARZCST,zlistevert,ikdebzint,XSPVAL)
1126                   ELSE IF (INDEX(YTYPEOUT(1:4),'A')/=0 .OR. INDEX(YTYPEOUT(1:4),'a')/=0) THEN
1127                     IF (.NOT. ALLOCATED(zlistevert3D)) THEN
1128                       ALLOCATE(zlistevert3D(SIZE(XZS,1),SIZE(XZS,2),SIZE(zlistevert)))
1129                       DO IZLIST=1,SIZE(zlistevert)
1130                         zlistevert3D(:,:,IZLIST)=XZS(:,:)+zlistevert(IZLIST)
1131                       ENDDO
1132                     ENDIF     
1133                     CALL SINTER(ZWORK3D,XZZ,ZVARZCST,zlistevert3D,ikdebzint,XSPVAL)
1134                   ELSE IF (INDEX(YTYPEOUT(1:4),'P')/=0 .OR. INDEX(YTYPEOUT(1:4),'p')/=0) THEN
1135                     CALL PINTER(ZWORK3D,0,XSPVAL,zlistevert,ZVARZCST,ZPABS)
1136                   ELSE IF (INDEX(YTYPEOUT(1:4),'H')/=0 .OR. INDEX(YTYPEOUT(1:4),'h')/=0) THEN
1137                     ZVARZCST(:,:,:)=ZWORK3D(:,:,:)
1138                   ELSE
1139                     print*,'** ERROR in vertical interpolations with ',YTYPEOUT
1140                   ENDIF
1141                   XVAR(:,:,:,J4,J5,J6)=ZVARZCST
1142                 END DO
1143               END DO
1144             END DO
1145             DEALLOCATE(ZVARSAVE,ZVARZCST,ZWORK3D)
1146             zmini=MINVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
1147             zmaxi=MAXVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
1148             print * ,' After vertical interpolation, min,max of the variable ',TRIM(CGROUP),'=', zmini,zmaxi
1149             ivarkdeb=1
1150             ivarkfin=inbvertz
1151             IF (ilocverbia >= 5 ) then
1152               print*,'ivarkdeb,ivarkfin= ',ivarkdeb,ivarkfin 
1153             ENDIF
1154           ENDIF
1155         ENDIF
1156         ! c. interpolation eventuelle sur l horizontale
1157         IF ( YOUTGRID(1:4) == 'LALO' ) THEN
1158           if (ilocverbia >= 0 ) then
1159             print *,'Translate to a regular lat lon grid '
1160           end if
1161           IF ( .NOT. ALLOCATED (ZNEWX) ) THEN
1162             IF ( IFLAGzcst == 1 ) THEN
1163               print*,'** no processing with VLEV at the first group'
1164               GOTO 99
1165             ELSE
1166             ! c.1. creation de la grille réguliere en lat lon
1167               if (ilocverbia >= 2 ) then
1168                 print *,'grid creation, size of XLON: ',SIZE(XLON,1),SIZE(XLON,2) 
1169               end if
1170               ! calcul des coord X Y des points de la grille lat-lon reguliere
1171               ! determine le maximum d espacement en lat et lon sur le zoom
1172               ALLOCATE(ZDIFFLON(SIZE(XLON,1)-1,SIZE(XLON,2)))
1173               ALLOCATE(ZDIFFLAT(SIZE(XLAT,1),SIZE(XLAT,2)-1))
1174               
1175               DO ii=1,SIZE(XLON,1)-1
1176                 DO jj=1,SIZE(XLON,2)
1177                      ZDIFFLON(ii,jj)=XLON(ii+1,jj)-XLON(ii,jj)
1178                 END DO
1179               END DO
1180
1181               DO ii=1,SIZE(XLAT,1)
1182                 DO jj=1,SIZE(XLAT,2)-1
1183                      ZDIFFLAT(ii,jj)=XLAT(ii,jj+1)-XLAT(ii,jj)
1184                 END DO
1185               END DO
1186 !              ZDELTALON=NINT(maxval(ZDIFFLON)*1000.)
1187 !              ZDELTALAT=NINT(maxval(ZDIFFLAT)*1000.)
1188                            ZDELTALON=maxval(ZDIFFLON)
1189               ZDELTALAT=maxval(ZDIFFLAT)
1190               DEALLOCATE(ZDIFFLON)
1191               DEALLOCATE(ZDIFFLAT)
1192               if (ZDELTALON == 0 .OR. ZDELTALAT == 0 ) THEN
1193                 print *,' error during ZDELTALON,ZDELTALAT computation=', ZDELTALON,ZDELTALAT
1194                 print *,'XLON(ivarideb+1,ivarjdeb)-XLON(ivarideb,ivarjdeb)'&
1195                         ,'XLON(ivarifin,ivarjfin)-XLON(ivarifin-1,ivarjfin)'&
1196                         ,'XLAT(ivarideb,ivarjdeb+1)-XLAT(ivarideb,ivarjdeb)'&
1197                         ,'XLAT(ivarifin,ivarjfin)-XLAT(ivarifin,ivarjfin-1)'
1198                 print *,XLON(ivarideb+1,ivarjdeb)-XLON(ivarideb,ivarjdeb)&
1199                        ,XLON(ivarifin,ivarjfin)-XLON(ivarifin-1,ivarjfin)&
1200                        ,XLAT(ivarideb,ivarjdeb+1)-XLAT(ivarideb,ivarjdeb)&
1201                        ,XLAT(ivarifin,ivarjfin)-XLAT(ivarifin,ivarjfin-1)           
1202                 print *, 'ivarideb+1,ivarjdeb,ivarifin-1,ivarjfin',ivarideb+1,ivarjdeb,ivarifin-1,ivarjfin
1203                 print *,'Verify the fields LAT LON of the FM file'
1204                 ALLOCATE(ZX(SIZE(XLAT,1),SIZE(XLAT,2)),ZY(SIZE(XLAT,1),SIZE(XLAT,2)))
1205                 ZX(1:SIZE(XZZ,1),1) = XXX(1:SIZE(XZZ,1),IGRID)
1206                 ZX(:,2:SIZE(XZZ,2)) = SPREAD(ZX(:,1),2,SIZE(XZZ,2)-1)
1207                 ZY(1,1:SIZE(XZZ,2)) = XXY(1:SIZE(XZZ,2),IGRID)
1208                 ZY(2:SIZE(XZZ,1),:) = SPREAD(ZY(1,:),1,SIZE(XZZ,1)-1)
1209                 CALL SM_LATLON(XLATORI,XLONORI,ZX,ZY,XLAT,XLON)
1210                 ZDELTALON=max(XLON(ivarideb+1,ivarjdeb)-XLON(ivarideb,ivarjdeb)&
1211                            ,XLON(ivarifin,ivarjfin)-XLON(ivarifin-1,ivarjfin))
1212                 ZDELTALAT=max(XLAT(ivarideb,ivarjdeb+1)-XLAT(ivarideb,ivarjdeb)&
1213                            ,XLAT(ivarifin,ivarjfin)-XLAT(ivarifin,ivarjfin-1))
1214                 print *,' After Model Grid computation: ZDELTALON,ZDELTALAT=', ZDELTALON,ZDELTALAT
1215               endif
1216               IDIM1=(maxval(XLON)-minval(XLON))/ZDELTALON
1217               IDIM2=(maxval(XLAT)-minval(XLAT))/ZDELTALAT
1218               ALLOCATE (ZNEWLAT(IDIM1,IDIM2),ZNEWLON(IDIM1,IDIM2) )
1219               if (ilocverbia >= 1 ) then
1220                 print*,' ZDELTALON,ZDELTALAT= ',ZDELTALON,ZDELTALAT
1221               endif
1222               if (ilocverbia >= 2 ) then
1223                 print*,' IDIM1,IDIM2= ',IDIM1,IDIM2
1224               endif
1225               ! depart de la nouvelle grille : coin Sud Ouest
1226               DO JI=1,IDIM1
1227                 ZNEWLON(JI,:)=minval(XLON) + (JI-1) *ZDELTALON
1228               ENDDO
1229               DO JJ=1,IDIM2
1230                 ZNEWLAT(:,JJ)=minval(XLAT) + (JJ-1) *ZDELTALAT
1231               ENDDO
1232               if (ilocverbia >= 4 ) then
1233                 print*, 'new lat lon grid=',ZNEWLAT(1,:)
1234                 print*, ZNEWLON(:,1)
1235               endif
1236
1237               ALLOCATE (ZNEWX(IDIM1,IDIM2))
1238               ALLOCATE (ZNEWY(IDIM1,IDIM2))
1239               CALL SM_XYHAT(XLATORI,XLONORI,ZNEWLAT,ZNEWLON,ZNEWX,ZNEWY)
1240               if (ilocverbia >= 4 ) then
1241                 ! XXX= XXHAT et XXY=XYHAT pour les 7 grilles
1242                 print*,' After SM_XYHAT old limits X: ', &
1243                        XXX(1,IGRID),XXX(SIZE(XVAR,1),IGRID)
1244                 print*,'                new limits X: ', &
1245                        ZNEWX(1,1),ZNEWX(IDIM1,IDIM2)
1246                 print*,'                old limits Y: ', &
1247                        XXY(1,IGRID),XXY(SIZE(XVAR,2),IGRID)
1248                 print*,'                new limits Y: ', &
1249                        ZNEWY(1,1),ZNEWY(IDIM1,IDIM2)
1250               endif
1251               if (ilocverbia >= 5 ) then
1252                 DO JI= 1,SIZE(XVAR,1) 
1253                   print*,'XXHAT ZNEWX',XXX(JI,IGRID),ZNEWX(JI,1),ZNEWX(JI,IDIM2)
1254                 ENDDO
1255                 DO JJ= 1,SIZE(XVAR,2) 
1256                   print*,'XYHAT ZNEWY',XXY(JJ,IGRID),ZNEWY(1,JJ),ZNEWX(IDIM1,JJ)
1257                 ENDDO
1258               endif
1259               ! calcul de la section de tableau correspondant au zoom
1260
1261 !===================================================================================================================
1262               I1=(maxval(XLON(ivarideb:ivarifin,ivarjdeb:ivarjfin)) &
1263                  -minval(XLON(ivarideb:ivarifin,ivarjdeb:ivarjfin)) )/ZDELTALON
1264               I2=(maxval(XLAT(ivarideb:ivarifin,ivarjdeb:ivarjfin)) &
1265                  -minval(XLAT(ivarideb:ivarifin,ivarjdeb:ivarjfin)) )/ZDELTALAT
1266               IZOOMIDEB=MAX(MIN(COUNT(ZNEWLON(:,1)<minval(XLON(ivarideb:ivarifin,ivarjdeb:ivarjfin))),IDIM1),1)
1267               IZOOMJDEB=MAX(MIN(COUNT(ZNEWLAT(1,:)<minval(XLAT(ivarideb:ivarifin,ivarjdeb:ivarjfin))),IDIM2),1)
1268               IZOOMIFIN=MAX(MIN(COUNT(ZNEWLON(:,1)<maxval(XLON(ivarideb:ivarifin,ivarjdeb:ivarjfin))),IDIM1),1)
1269               IZOOMJFIN=MAX(MIN(COUNT(ZNEWLAT(1,:)<maxval(XLAT(ivarideb:ivarifin,ivarjdeb:ivarjfin))),IDIM2),1)
1270 !=====================================================================================================================
1271
1272
1273               if (ilocverbia >= 2 ) then
1274                 print*,' ZOOM along i in the LON-LAT grid: ', &
1275                        IZOOMIDEB,IZOOMIFIN,I1
1276                 print*,'            j                    : ', &
1277                        IZOOMJDEB,IZOOMJFIN,I2
1278               endif
1279             ENDIF
1280           ENDIF ! fin grille ZNEWX deja allouee
1281           ! c.2. interpolation sur la nouvelle grille
1282           IF( IFLAGzcst/= 1 .AND. (NREADIH-NREADIL)>0 .AND. (NREADJH-NREADJL)>0 )THEN
1283             ! interpolation vers la nouvelle grille réguliere en lat lon 
1284             !sauf la grille verticale definie en niveaux Z et champs 1D
1285             if (ilocverbia >= 1 ) then
1286               print*,' interpolation for the variable  ',trim(CGROUP)
1287             end if
1288             allocate(ZWORK3D(SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3)))
1289             allocate(ZWORK3D2(IDIM1,IDIM2,SIZE(XVAR,3)))
1290             ! stockage des champs interpoles dans la nouvelle grille
1291             if (allocated (ZVARSAVE)) DEALLOCATE(ZVARSAVE)
1292             allocate(ZVARSAVE(IDIM1,IDIM2,SIZE(XVAR,3),&
1293                      SIZE(XVAR,4),SIZE(XVAR,5),SIZE(XVAR,6)))
1294             ! boucle sur les dimensions 4 5 6
1295             DO J6=ivarprocinf,ivarprocsup
1296               DO J5=ivartrajinf,ivartrajsup
1297                 DO J4=ivartinf,ivartsup
1298                   ZWORK3D(:,:,:)=XVAR(:,:,:,J4,J5,J6)
1299                   if (ilocverbia >= 2 ) then
1300                     print *,'before HOR_INTERP_4PTS J4,J5,J6=', J4,J5,J6
1301                   end if
1302                   CALL HOR_INTERP_4PTS(XXX(:,IGRID),XXY(:,IGRID),ZWORK3D, &
1303                                        ZNEWX,ZNEWY,ZWORK3D2)
1304                   ZVARSAVE(:,:,:,J4,J5,J6)=ZWORK3D2(:,:,:)
1305                 END DO
1306               END DO
1307             END DO
1308             ! resultat dans XVAR passe en module
1309             DEALLOCATE (XVAR)
1310             ALLOCATE(XVAR(IDIM1,IDIM2,SIZE(ZVARSAVE,3),&
1311                     SIZE(ZVARSAVE,4),SIZE(ZVARSAVE,5),SIZE(ZVARSAVE,6)))
1312             XVAR=XSPVAL
1313             XVAR(:,:,:,ivartinf:ivartsup,ivartrajinf:ivartrajsup,ivarprocinf:ivarprocsup)= &
1314             ZVARSAVE(:,:,:,ivartinf:ivartsup,ivartrajinf:ivartrajsup,ivarprocinf:ivarprocsup)
1315             DEALLOCATE (ZVARSAVE)
1316             zmini=MINVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
1317             zmaxi=MAXVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
1318             print * ,' After horizontal interpolation, min,max of the variable ',TRIM(CGROUP),'=', zmini,zmaxi
1319             if (ilocverbia >= 2 ) then
1320               print*, 'After HOR_INTERP_4PTS all the dim 4,5,6'
1321             endif
1322             IF (allocated(ZWORK3D)) DEALLOCATE(ZWORK3D)
1323             IF (allocated(ZWORK3D2)) DEALLOCATE(ZWORK3D2)
1324           ENDIF
1325         ENDIF
1326         ! d. ecriture des donnees au format cdl ou llz/llp
1327         IF ( YOUTGRID(1:4) == 'LALO' ) THEN
1328           IF ( IFLAGzcst /= 1 ) THEN
1329             ivarideb=IZOOMIDEB
1330             ivarifin=IZOOMIFIN
1331             ivarjdeb=IZOOMJDEB
1332             ivarjfin=IZOOMJFIN
1333           ENDIF
1334           SELECT CASE(YTYPEOUT(1:4))
1335           CASE('LLZV','llzv','LLPV','llpv','LLAV','llav')
1336             IF (allocated(ZWORK3D)) DEALLOCATE(ZWORK3D)
1337             ALLOCATE(ZWORK3D(size(XVAR,1),size(XVAR,2),size(XVAR,3)))
1338             IF (SIZE(XVAR,3)==inbvertz) THEN
1339               ZWORK3D(1,1,:)=zlistevert
1340             ELSE
1341               ZWORK3D(1,1,:)=XSPVAL
1342             ENDIF
1343             CALL WRITELLHV(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1344                        ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1345                        ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1346                        CGROUP,YFILEIN,YFLAGWRITE,YTYPEOUT,&
1347                        ilocverbia,iret,PLON=ZNEWLON,PLAT=ZNEWLAT,&
1348                        PALT=ZWORK3D)       
1349             if (ilocverbia > 0 ) then
1350               print*,'WRITELLHV LALO return= ', YTYPEOUT,'= ',iret
1351             end if
1352             DEALLOCATE(ZWORK3D)
1353         !
1354           CASE('KCDL','ZCDL','PCDL')
1355            YGROUP=ADJUSTL(ADJUSTR(CGROUP)//ADJUSTL(YK))
1356            CALL WRITECDL(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1357                         ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1358                         ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup, &
1359                         YGROUP,YFILEIN,YFLAGWRITE,YOUTGRID,YSUFFIX_file, &
1360                         ilocverbia,iret,PGRIDX=ZNEWLON(:,1),PGRIDY=ZNEWLAT(1,:))
1361            IF (ilocverbia >= 1 ) print *,' counter of added fields=',inetadd
1362            if ( inetadd == 0) then
1363             IF( SIZE(XZZ,3)<=1 .OR. SIZE(XZZ,2)<=1 .OR. SIZE(XZZ,1)<=1 ) THEN
1364             ! VLEV, LON, LAT et chps 2D ne passent pas cette partie 
1365              print *,' *****The program could not add the VLEV 3Dfield to the netcdf file****'
1366             ELSE
1367              print *,' The program adds the VLEV 3Dfield to the netcdf file'
1368              YGROUP_OLD=CGROUP(1:13)
1369              CGROUP='VLEV'
1370              inetadd=inetadd+1
1371              YFLAGWRITE='OLD'
1372              ino_init_zoom=1
1373              GO TO 77
1374             ENDIF
1375            endif
1376            if ( inetadd == 1 .AND. YOUTGRID(1:4) == 'CONF' )THEN
1377              print *,' The program adds the LAT 3Dfield to the netcdf file'      
1378              CGROUP='LAT'
1379              inetadd=inetadd+1
1380              YFLAGWRITE='OLD'
1381              ino_init_zoom=1
1382              GO TO 77
1383            endif
1384            if ( inetadd == 2 .AND. YOUTGRID(1:4) == 'CONF' )THEN
1385              print *,' The program adds the LON 3Dfield to the netcdf file'      
1386              CGROUP='LON'
1387              inetadd=inetadd+1
1388              YFLAGWRITE='OLD'
1389              ino_init_zoom=1
1390              GO TO 77
1391            endif
1392            
1393           END SELECT
1394         ELSE ! pas d interpolation horizontale
1395           SELECT CASE(YTYPEOUT(1:4))
1396           CASE('LLZV','llzv','LLPV','llpv','IJZV','jizv','IJPV','jipv','LLAV','llav')
1397             IF (SIZE(XVAR,3)==inbvertz) THEN  ! champ 3D
1398               IF (allocated(ZWORK3D)) DEALLOCATE(ZWORK3D)
1399               ALLOCATE(ZWORK3D(size(XVAR,1),size(XVAR,2),size(XVAR,3)))
1400               ZWORK3D(1,1,:)=zlistevert
1401               CALL WRITELLHV(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1402                        ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1403                        ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1404                        CGROUP,YFILEIN,YFLAGWRITE,YTYPEOUT,&
1405                        ilocverbia,iret,&
1406                        PALT=ZWORK3D)       
1407             ELSE                              ! champ 2D
1408               IF((YTYPEOUT(3:3)=='z').OR.(YTYPEOUT(3:3)=='p')) YTYPEOUT3='h'
1409               IF((YTYPEOUT(3:3)=='Z').OR.(YTYPEOUT(3:3)=='P')) YTYPEOUT3='H'
1410               IF((YTYPEOUT(3:3)=='a').OR.(YTYPEOUT(3:3)=='A')) YTYPEOUT3='H'
1411               CALL WRITELLHV(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1412                        ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1413                        ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup, &
1414                        CGROUP,YFILEIN,YFLAGWRITE, &
1415                        YTYPEOUT(1:2)//YTYPEOUT3//YTYPEOUT(4:4), &
1416                        ilocverbia,iret)
1417             ENDIF
1418             if (ilocverbia > 0 ) then
1419               print*,' WRITELLHV for ', YTYPEOUT,', return value= ',iret
1420             end if
1421             IF (allocated(ZWORK3D)) DEALLOCATE(ZWORK3D)
1422         !
1423           CASE('KCDL','ZCDL','PCDL')
1424             YGROUP=ADJUSTL(ADJUSTR(CGROUP)//ADJUSTL(YK))
1425             CALL WRITECDL(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1426                         ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1427                         ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup, &
1428                         YGROUP,YFILEIN,YFLAGWRITE,YOUTGRID,YSUFFIX_file, &
1429                         ilocverbia,iret,PGRIDX=XXX(:,IGRID),PGRIDY=XXY(:,IGRID))
1430             IF (ilocverbia >= 1 ) print *,' counter of added fields=',inetadd
1431            if ( inetadd == 0) then
1432              if (ivarkdeb == ivarkfin .AND. ivarkdeb == 1 ) THEN
1433                print *, 'No VLEV field for only one vertical position'
1434              else
1435              print *,' The program adds the VLEV 3Dfield to the netcdf file'
1436              YGROUP_OLD=CGROUP(1:13)
1437              CGROUP='VLEV'
1438              inetadd=inetadd+1
1439              YFLAGWRITE='OLD'
1440              ino_init_zoom=1
1441              GO TO 77
1442              endif
1443            endif
1444            if ( inetadd == 1 .AND. YOUTGRID(1:4) == 'CONF' )THEN
1445              if (ivarideb /= ivarifin .AND. ivarjdeb /= ivarjfin ) THEN
1446
1447              print *,' The program adds the LAT 3Dfield to the netcdf file'      
1448              CGROUP='LAT'
1449              inetadd=inetadd+1
1450              ino_init_zoom=1
1451              GO TO 77
1452              else
1453                print *, ' No LAT field for only one location', ivarideb,ivarifin,ivarjdeb,ivarjfin
1454              endif
1455            endif
1456            if ( inetadd == 2 .AND. YOUTGRID(1:4) == 'CONF' )THEN
1457              if (ivarideb /= ivarifin .AND. ivarjdeb /= ivarjfin ) THEN      
1458              print *,' The program adds the LON 3Dfield to the netcdf file'      
1459              CGROUP='LON'
1460              inetadd=inetadd+1
1461              ino_init_zoom=1
1462              GO TO 77
1463              else
1464                print *, ' No LON field for only one location', ivarideb,ivarifin,ivarjdeb,ivarjfin
1465              endif
1466            endif
1467           END SELECT              
1468         ENDIF
1469         ! retour a XZZ pour NGRID a 4 (cf readvar)
1470         CALL COMPCOORD_FORDIACHRO(4)
1471 !============================================   
1472    CASE('ZGRB','PGRB','AGRB')
1473         IF(SIZE(XVAR,3)==1) THEN
1474            LVAR2D=.TRUE.
1475         ENDIF
1476         ! replace field at mass points
1477         IF ( CGROUP /= 'VLEV' ) THEN
1478           If (ALLOCATED(ZWORK3D))DEALLOCATE(ZWORK3D)
1479           If (ALLOCATED(ZWORK3D2))DEALLOCATE(ZWORK3D2)
1480           ALLOCATE(ZWORK3D(SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3)))
1481           ALLOCATE(ZWORK3D2(SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3)))
1482           DO J6=ivarprocinf,ivarprocsup
1483             IGRID=NGRIDIA(J6)
1484             IF(SIZE(XVAR,3)/=1 .OR. IGRID/=4) THEN 
1485               ! pas d interpolation verticale pour champ 2D
1486               DO J5=ivartrajinf,ivartrajsup
1487                 DO J4=ivartinf,ivartsup
1488                   ZWORK3D(:,:,:)=XVAR(:,:,:,J4,J5,J6)
1489                   print'(A29,3(X,I4))',' mass point grid for J4,J5,J6=',J4,J5,J6
1490                   CALL CHANGE_A_GRID(ZWORK3D,IGRID,ZWORK3D2)
1491                   ! IGRID=1 en sortie de change_a_grid
1492                   NGRIDIA(J6)=IGRID
1493                   XVAR(:,:,:,J4,J5,J6)=ZWORK3D2(:,:,:)               
1494                 ENDDO
1495               ENDDO
1496             ENDIF
1497           ENDDO
1498           DEALLOCATE(ZWORK3D,ZWORK3D2)
1499         ENDIF
1500         !
1501         ! a. reinit avant ecriture de la grille verticale correspondant a la
1502         !grille de masse sur laquelle le champ a ete interpole
1503         IFLAGzcst=0
1504         IF (IND_VERT/=0) THEN
1505           ! b. interpolation eventuelle selon la verticale 
1506           IF( SIZE(XVAR,3)>1 .AND. SIZE(XVAR,2)>1 .AND. SIZE(XVAR,1)>1 ) THEN
1507             ! VLEV, LON, LAT et chps 2D ne passent pas cette partie 
1508             if (ilocverbia >= 0 ) then
1509               print*,' Interpolations on ',inbvertz,' ', &
1510                      YTYPEOUT(IND_VERT:IND_VERT),'-levels'
1511             endif
1512             if (ilocverbia >= 1 .AND. IND_VERT/=0) THEN
1513               print*,'levels= ',zlistevert 
1514             endif
1515             ALLOCATE(ZVARSAVE(size(XVAR,1),size(XVAR,2),size(XVAR,3),   &
1516                               size(XVAR,4),size(XVAR,5),size(XVAR,6))   )
1517             ZVARSAVE=XVAR
1518             ALLOCATE(ZWORK3D(SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3)))
1519             ALLOCATE(ZVARZCST(SIZE(XVAR,1),SIZE(XVAR,2),inbvertz))
1520             DEALLOCATE(XVAR)
1521             ALLOCATE(XVAR(SIZE(ZVARSAVE,1),SIZE(ZVARSAVE,2),SIZE(ZVARZCST,3),&
1522                           size(ZVARSAVE,4),size(ZVARSAVE,5),size(ZVARSAVE,6)))
1523             DO J6=ivarprocinf,ivarprocsup
1524               IGRID=NGRIDIA(J6)
1525               ! init du tableau des altitudes  XZZ pour la grille= IGRID
1526               CALL COMPCOORD_FORDIACHRO(IGRID)
1527               DO J5=ivartrajinf,ivartrajsup
1528                 DO J4=ivartinf,ivartsup
1529                   ZWORK3D(:,:,:)=ZVARSAVE(:,:,:,J4,J5,J6)
1530                   ikdebzint=2
1531                   IF (INDEX(YTYPEOUT(1:4),'Z')/=0 .OR. INDEX(YTYPEOUT(1:4),'z')/=0) THEN
1532                     CALL ZINTER(ZWORK3D,XZZ,ZVARZCST,zlistevert,ikdebzint,XSPVAL)
1533                   ELSE IF (INDEX(YTYPEOUT(1:4),'A')/=0 .OR. INDEX(YTYPEOUT(1:4),'a')/=0) THEN
1534                     IF (.NOT. ALLOCATED(zlistevert3D)) THEN
1535                       ALLOCATE(zlistevert3D(SIZE(XZS,1),SIZE(XZS,2),SIZE(zlistevert)))
1536                       DO IZLIST=1,SIZE(zlistevert)
1537                         zlistevert3D(:,:,IZLIST)=XZS(:,:)+zlistevert(IZLIST)
1538                       ENDDO
1539                     ENDIF     
1540                     CALL SINTER(ZWORK3D,XZZ,ZVARZCST,zlistevert3D,ikdebzint,XSPVAL)
1541                   ELSE IF (INDEX(YTYPEOUT(1:4),'P')/=0 .OR. INDEX(YTYPEOUT(1:4),'p')/=0) THEN
1542                     CALL PINTER(ZWORK3D,0,XSPVAL,zlistevert,ZVARZCST,ZPABS)
1543                   ELSE IF (INDEX(YTYPEOUT(1:4),'H')/=0 .OR. INDEX(YTYPEOUT(1:4),'h')/=0) THEN
1544                     ZVARZCST(:,:,:)=ZWORK3D(:,:,:)
1545                   ELSE
1546                     print*,'** ERROR in vertical interpolations with ',YTYPEOUT
1547                   ENDIF
1548                   XVAR(:,:,:,J4,J5,J6)=ZVARZCST
1549                 END DO
1550               END DO
1551             END DO
1552             DEALLOCATE(ZVARSAVE,ZVARZCST,ZWORK3D)
1553             zmini=MINVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
1554             zmaxi=MAXVAL(XVAR(:,:,:,:,:,:),MASK=XVAR(:,:,:,:,:,:)/=XSPVAL)
1555             print * ,' After vertical interpolation, min,max of the variable ',TRIM(CGROUP),'=', zmini,zmaxi
1556             ivarkdeb=1
1557             ivarkfin=inbvertz
1558             IF (ilocverbia >= 5 ) then
1559               print*,'ivarkdeb,ivarkfin= ',ivarkdeb,ivarkfin 
1560             ENDIF
1561           ENDIF
1562         ENDIF
1563         ! c. interpolation eventuelle sur l horizontale
1564         IF ( YOUTGRID(1:4) == 'LALO' ) THEN
1565           ZLATLON(1)=MAXVAL(XLAT)*1000.
1566           ZLATLON(2)=MINVAL(XLAT)*1000.
1567           ZLATLON(3)=MINVAL(XLON)*1000.
1568           ZLATLON(4)=MAXVAL(XLON)*1000.
1569       
1570           IF (ZJFIN /=0 .AND. ZJFIN/=-1) ZLATLON(1)=zjfin*1000.
1571           IF (ZJDEB /=0 .AND. ZJDEB/=-1) ZLATLON(2)=zjdeb*1000.
1572           IF (ZIDEB /=0 .AND. ZIDEB/=-1) ZLATLON(3)=zideb*1000.
1573           IF (ZIFIN /=0 .AND. ZIFIN/=-1) ZLATLON(4)=zifin*1000.
1574          
1575           ILATLON(:)=ZLATLON(:)
1576           
1577           IF (ILATLON(1)>  ZLATLON(1)) ILATLON(1)=ILATLON(1)-1
1578           IF (ILATLON(2)<  ZLATLON(2)) ILATLON(2)=ILATLON(2)+1
1579           IF (ILATLON(3)<  ZLATLON(3)) ILATLON(3)=ILATLON(3)+1
1580           IF (ILATLON(4)>  ZLATLON(4)) ILATLON(4)=ILATLON(4)-1
1581  
1582           ZLATLON=ILATLON
1583           IF (CFIXRESOL=="y") THEN
1584             INX=(ZLATLON(4)-ZLATLON(3))/ZDX_GRB +1 
1585             INY=(ZLATLON(1)-ZLATLON(2))/ZDY_GRB +1 
1586             ZCONTROL=ZLATLON(3)+(INX-1)*ZDX_GRB
1587             IF (ZCONTROL/=ZLATLON(4)) THEN
1588               print*,"warning : need to change E longitude"
1589               ZLATLON(4)=ZCONTROL
1590               print*,"lon min" ,ZLATLON(3)
1591               print*,"lon max" ,ZLATLON(4)
1592               print*,"dx",ZDX_GRB
1593             ENDIF
1594             ZCONTROL=ZLATLON(2)+(INY-1)*ZDY_GRB
1595             IF (ZCONTROL/=ZLATLON(1)) THEN
1596               print*,"warning : need to change N latitude"
1597               ZLATLON(1)=ZCONTROL
1598               print*,"lat min" ,ZLATLON(2)
1599               print*,"lat max" ,ZLATLON(1)
1600               print*,"dy",ZDY_GRB
1601             ENDIF
1602             print*,INX,INY,ZLATLON,ZDX_GRB,ZDY_GRB
1603           ELSE
1604             CALL INI2LALO(ZLATLON,INX,INY)
1605           ENDIF
1606           ALLOCATE(ZVARSAVE(size(XVAR,1),size(XVAR,2),size(XVAR,3),   &
1607                               size(XVAR,4),size(XVAR,5),size(XVAR,6))   )
1608           ZVARSAVE=XVAR
1609           ALLOCATE(ZWORK3D(SIZE(XVAR,1),SIZE(XVAR,2),SIZE(XVAR,3)))
1610           ALLOCATE(ZVARZCST(INX,INY,size(XVAR,3)))
1611           DEALLOCATE(XVAR)
1612           ALLOCATE(XVAR(INX,INY,SIZE(ZVARSAVE,3),&
1613                           size(ZVARSAVE,4),size(ZVARSAVE,5),size(ZVARSAVE,6)))
1614
1615           DO J6=ivarprocinf,ivarprocsup
1616             DO J5=ivartrajinf,ivartrajsup
1617               DO J4=ivartinf,ivartsup
1618                 ZWORK3D(:,:,:)=ZVARSAVE(:,:,:,J4,J5,J6)
1619                 CALL INT2LALO('BILI',ZWORK3D,ZLATLON,XSPVAL,ZVARZCST)
1620                 XVAR(:,:,:,J4,J5,J6)=ZVARZCST
1621               END DO
1622             END DO
1623           END DO
1624           DEALLOCATE(ZVARSAVE,ZVARZCST,ZWORK3D)
1625         ENDIF
1626        print*," ZLATLON apres INT2lalo",ZLATLON
1627                  ! d. ecriture des donnees au format  GRIB
1628         IF ( YOUTGRID(1:4) == 'LALO' ) THEN
1629           IF ( IFLAGzcst /= 1 ) THEN
1630             ivarideb=1
1631             ivarifin=SIZE(XVAR,1)
1632             ivarjdeb=1
1633             ivarjfin=SIZE(XVAR,2)
1634           ENDIF
1635           IF (LLEVEL2D) THEN
1636                  CALL WRITEGRIB(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1637                        ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1638                        ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1639                        CGROUP,YFILEIN,YFLAGWRITE,YOUTGRID,YTYPEOUT,&
1640                        ilocverbia,iret,ICODCOD,&
1641                        zlistevert,LVAR2D,KLEVEL2D=ILEVEL2D,PLATLON=ZLATLON)        
1642           ELSE
1643                  CALL WRITEGRIB(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1644                        ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1645                        ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1646                        CGROUP,YFILEIN,YFLAGWRITE,YOUTGRID,YTYPEOUT,&
1647                        ilocverbia,iret,ICODCOD,&
1648                        zlistevert,LVAR2D,PLATLON=ZLATLON)        
1649           ENDIF 
1650           if (ilocverbia > 0 ) then
1651              print*,'WRITEGRIB LALO return= ', YTYPEOUT,'= ',iret
1652           end if
1653         ELSE ! pas d interpolation horizontale (CONF)
1654            IF (LCARTESIAN) THEN
1655                  PRINT*,"===================================="
1656                  PRINT*,"WARNING : WITH LCARTESIAN=TRUE PLEASE ASK LALO"
1657                  PRINT*,"===================================="
1658                  STOP      
1659            ENDIF
1660            IF (LLEVEL2D) THEN
1661                  CALL WRITEGRIB(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1662                        ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1663                        ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1664                        CGROUP,YFILEIN,YFLAGWRITE,YOUTGRID,YTYPEOUT,&
1665                        ilocverbia,iret,ICODCOD,zlistevert,LVAR2D,KLEVEL2D=ILEVEL2D)
1666            ELSE
1667                  CALL WRITEGRIB(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1668                        ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1669                        ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1670                        CGROUP,YFILEIN,YFLAGWRITE,YOUTGRID,YTYPEOUT,&
1671                        ilocverbia,iret,ICODCOD,zlistevert,LVAR2D)      
1672            ENDIF         
1673            if (ilocverbia > 0 ) then
1674                  print*,'WRITEGRIB CONF return= ', YTYPEOUT,'= ',iret
1675            end if
1676         ENDIF
1677         ! retour a XZZ pour NGRID a 4 (cf readvar)
1678         CALL COMPCOORD_FORDIACHRO(4)
1679      END SELECT
1680      ! indiquera aux routines d ecriture que le fichier courant est deja ouvert
1681      YFLAGWRITE='OLD'
1682   ! 
1683   ELSE   ! iret /=0
1684     print *, ' READVAR return= ',iret
1685   ENDIF  
1686 END DO ! boucle champ a traiter
1687 !
1688 !
1689 !---------------------------------------------------------------------------
1690 !
1691 !*       4.    CLOSURE OF OUTPUT FILE
1692 !              ----------------------
1693 !
1694 !pour clore le traitement meme si la liste des champs est non terminee par END
1695 88 CONTINUE
1696 !
1697 IF (ALLOCATED(ZNEWX))   DEALLOCATE(ZNEWX,ZNEWY)
1698 IF (ALLOCATED(ZNEWLAT)) DEALLOCATE(ZNEWLAT,ZNEWLON)
1699 IF (ALLOCATED(ZWORK2D)) DEALLOCATE(ZWORK2D,ZWORK2D2)
1700 !
1701 PRINT*, 'END ->  Close the output file'
1702 YFLAGWRITE='CLO'
1703 SELECT CASE(YTYPEOUT(1:4))
1704   CASE('DIAC')
1705     CALL WRITEVAR(ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
1706                   ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,  &
1707                  CGROUP,YFILEIN,YFLAGWRITE,'2  ',ilocverbia,iret)
1708   CASE('LLHV','llhv','LLZV','llzv','LLPV','llpv','LLAV','llav',&
1709           'IJHV','IJZV','IJPV','jihv','jizv','jipv')             
1710     CALL WRITELLHV(ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
1711                  ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1712                  CGROUP,YFILEIN,YFLAGWRITE,YTYPEOUT,ilocverbia,iret)      
1713   CASE('KCDL','ZCDL','PCDL')
1714     CALL WRITECDL(ivarideb,ivarifin,ivarjdeb,ivarjfin,ivarkdeb,ivarkfin,&
1715                   ivartinf,ivartsup,ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1716                   CGROUP,YFILEIN,YFLAGWRITE,YOUTGRID,YSUFFIX_file,      &
1717                   ilocverbia,iret,PGRIDX=XXX(:,IGRID),PGRIDY=XXY(:,IGRID))
1718   CASE('ZGRB','PGRB','AGRB')
1719     CALL WRITEGRIB(ivarideb,ivarifin,ivarjdeb,ivarjfin, &
1720                    ivarkdeb,ivarkfin,ivartinf,ivartsup, &
1721                    ivartrajinf,ivartrajsup,ivarprocinf,ivarprocsup,&
1722                    CGROUP,YFILEIN,YFLAGWRITE,YOUTGRID,YTYPEOUT,&
1723                    ilocverbia,iret,ICODCOD,zlistevert,LVAR2D)   
1724   CASE DEFAULT
1725     PRINT*, 'Closure of output type ',YTYPEOUT ,' not coded'
1726 END SELECT
1727 !
1728 !-------------------------------------------------------------------------------
1729 !
1730 !*       5.    END
1731 !              ---
1732 !
1733 99 CONTINUE
1734 PRINT*, 'Delete the links if necessary'
1735 YDUMMYFILE=''
1736 CALL CREATLINK(' ',YDUMMYFILE,'CLEAN',ILOCVERBIA)
1737 PRINT*, 'The file ',TRIM(YLUDIR),' stores all the input directives '
1738 PRINT*, ' you must give a new name to use it again'
1739 CLOSE(ILUDIR)
1740 !
1741 !-------------------------------------------------------------------------------
1742 !
1743 CONTAINS 
1744 !
1745 !------------------------------------------------------------------------------
1746 !
1747 SUBROUTINE TAB2SPACE(HTEXT)
1748 IMPLICIT NONE
1749 CHARACTER(len=*),INTENT(INOUT) :: HTEXT
1750
1751 CHARACTER, PARAMETER :: YPTAB = CHAR(9) ! TAB character is ASCII : 9
1752 CHARACTER, PARAMETER :: YPCOM = CHAR(44)! COMma character is ASCII : 44
1753 INTEGER              :: JI
1754
1755 DO JI=1,LEN_TRIM(HTEXT)
1756   IF (HTEXT(JI:JI)==YPTAB .OR. HTEXT(JI:JI)==YPCOM) HTEXT(JI:JI) = ' '
1757 END DO
1758 END SUBROUTINE TAB2SPACE
1759 !------------------------------------------------------------------------------
1760
1761 END PROGRAM EXTRACTDIA
1762 !