Philippe 23/02/2016: lfi2cdf: modif for JPHEXT/=1
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / DIAPRO / compcoord_fordiachro.f90
1 !     ######spl
2       SUBROUTINE COMPCOORD_FORDIACHRO(KGRID)
3 !     ######################################
4 !
5 !!****  *COMPCOORD_FORDIACHRO* - Computes gridpoint locations,
6 !!                     meshsizes and topography 
7 !!                    for all the possible grids, and true altitude where 
8 !!                    required.
9 !!
10 !!    PURPOSE
11 !!    -------
12 !       When called for the first time (KGRID=0), COMPCOORD_FORDIACHRO returns for 
13 !     the 7 possible grid locations:
14 !        - XHAT, YHAT, ZHAT values (meters) stored in:
15 !              XXX(:,1:7), XXY(:,1:7), XXZ(:,1:7)
16 !        - meshsizes values (meters):
17 !              XXDXHAT(:,1:7), XXDYHAT(:,1:7)
18 !        - topography altitudes values (meters):
19 !              XXZS(:,:,1:7)
20 !
21 !       When called subsequently (0<KGRID<8), COMPCOORD_FORDIACHRO returns the true
22 !     gridpoint altitude (meters) corresponding to the requested KGRID value 
23 !     in the XZZ(:,:,:) array.
24 !
25 !!**  METHOD
26 !!    ------
27 !!      Temporary arrays are allocated to store the grid point characteristics
28 !!    and de-allocated on exit. The 3D gridpoints locations are linearly
29 !!    interpolated to the expected grid location from their respective
30 !!    nominal locations. Altitudes are interpolated for the w-grid values,
31 !!    which are obtained directly from the Gal-Chen Sommerville formula.
32 !!      For XXX, XXY, XXZ, XXDXHAT, XXDYHAT, XXZS the last index is the grid
33 !!    selector KGRID ranging from 1 to 7 as follows:
34 !!    1 -> Mass grid, 
35 !!    2 -> U grid, 
36 !!    3 -> V grid, 
37 !!    4 -> W grid, 
38 !!    5 -> Vertical vorticity grid, 
39 !!    6 -> y-component vorticity grid,
40 !!    7 -> x-component vorticity grid.
41 !!    all the 7 values are prepared one for all in this subroutine and passed
42 !!    to the general TRACE environment to be used in the display process.
43 !!
44 !!      For the XZZ array the last index is the z direction one, not the grid 
45 !!    selector one.
46 !! 
47 !!
48 !!    EXTERNAL
49 !!    --------
50 !!      None
51 !!
52 !!    IMPLICIT ARGUMENTS
53 !!    ------------------
54 !!      Module MODD_COORD      : declares gridpoint coordinates (TRACE use)
55 !!
56 !!       XXX,XXY,XXZ   : coordinate values for all the MESO-NH grids
57 !!       XXDXHAT,XXDYHAT,XXDZHAT:   meshsize values for all the MESO-NH grids
58 !!       XXZS          : topography values for all the MESO_NH grids
59 !!
60 !!      Module MODD_GRID1      : declares grid variables (Model module)
61 !!
62 !!       XXHAT, XYHAT : x, y in the conformal or cartesian plane
63 !!       XZHAT        : Gal-Chen z level
64 !!       XZS          : topography zs
65 !!       XZZ          : true gridpoint z altitude
66 !! 
67 !!      Module MODD_DIM1       : Contains dimensions
68 !!
69 !!         NIMAX,NJMAX,NKMAX :  x, y, and z array dimensions
70 !!
71 !!      Module MODD_PARAMETERS : Contains array border depths
72 !!
73 !!         JPHEXT   : Horizontal external points number
74 !!         JPVEXT   : Vertical external points number
75 !!
76 !!
77 !!    REFERENCE
78 !!    ---------
79 !!
80 !!      MESO-NH User's Manual, TRACE Post Processing sections, Version 1.0:
81 !!       + Book1: Concepts and Fundamentals, to appear in 1994;
82 !!       + Book2: Technical Reference and Flowcharts, to appear in 1994;
83 !!       + Book3: Tutorial, November 1994.
84 !!
85 !!     The 7 MESO-NH grid types are defined in:
86 !!      - Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
87 !!        commun CNRM-LA, specifications techniques",
88 !!        Note CNRM/GMME, 26, 139p, (pages 39 to 43).
89 !!
90 !!      - Fischer C., 1994, "File structure and content in the Meso-NH
91 !!        model", Meso-nh internal note, CNRM/GMME,  July 5.
92 !!
93 !!
94 !!    AUTHOR
95 !!    ------
96 !!      J. Duron    * Laboratoire d'Aerologie *
97 !!      
98 !!
99 !!    MODIFICATIONS
100 !!    -------------
101 !!      Original       06/06/94
102 !!      Updated   PM   02/12/94
103 !-------------------------------------------------------------------------------
104 !
105 !*       0.    DECLARATIONS
106 !              ------------
107 !
108 USE MODD_COORD
109 USE MODD_DIM1
110 USE MODD_CONF
111 USE MODD_GRID1
112 USE MODD_PARAMETERS
113 USE MODD_MEMCV
114 USE MODD_RESOLVCAR
115 !
116 USE MODI_VERT_COORD
117 !
118 IMPLICIT NONE
119 !
120 !*       0.1   Local variables declarations
121 !
122 INTEGER           :: IIU, IJU, IKU
123
124 INTEGER           :: IIB, IJB, IKB
125
126 INTEGER           :: IIE, IJE, IKE
127 !
128 ! Calcul des X, Y, Z aux points de masse
129 REAL,DIMENSION(:),ALLOCATABLE   :: ZXMASS, ZYMASS, ZZMASS 
130
131 REAL,DIMENSION(:),ALLOCATABLE   :: ZXTEM, ZYTEM, ZZTEM,  &
132                                    ZDXTEM, ZDYTEM, ZDZTEM
133
134 REAL,DIMENSION(:,:),ALLOCATABLE :: ZZSTEM
135
136 REAL,DIMENSION(:,:),ALLOCATABLE,SAVE   :: ZSCOEF
137
138 REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: ZSZ
139
140 REAL,SAVE :: ZSH
141
142 INTEGER :: JGRIDLOOP,  KGRID,  &
143            OKZXTEM, OKZYTEM, OKZZTEM, OKZDXTEM, OKZDYTEM, OKZDZTEM,  &
144            OKZZSTEM, OKZSCOEF, OKZXMASS, OKZYMASS, OKZZMASS, OKXXZS, OKXZZ, OKZSZ
145 !
146 !-------------------------------------------------------------------------------
147 !
148 !*       1.    ARRAY DIMENSIONS SETTOING
149 !              -------------------------
150 !
151 if(nverbia > 0)then
152 if (LEN(CDIRCUR) .LT. 500)THEN
153 print *,' **COMPCOORD KGRID DIRCUR ',KGRID,CDIRCUR(1:LEN_TRIM(CDIRCUR))
154 endif
155 endif
156 IIU=NIMAX+2*JPHEXT
157 IJU=NJMAX+2*JPHEXT
158 IKU=NKMAX+2*JPVEXT
159 IF(CSTORAGE_TYPE == 'PG' .OR. CSTORAGE_TYPE == 'SU')THEN
160   IKU=1
161 ENDIF
162 IIB=1+JPHEXT
163 IJB=1+JPHEXT
164 IKB=1+JPVEXT
165 IIE=IIU-JPHEXT
166 IJE=IJU-JPHEXT
167 IKE=IKU-JPVEXT
168 !
169 !-------------------------------------------------------------------------------
170 !
171 !*       2.   CALCULATIONS PERFORMED FOR THE FIRST CALL
172 !              ----------------------------------------
173 !
174 ! Test on KGRID selects processing mode:
175 ! . KGRID=0 --> X, Y, Z + meshsizes + topography computed for ALL the
176 !   possible grid geometry.
177 ! . 0<KGRID<8 --> true altitude computed for the KGRID gridpoints
178
179 IF(KGRID==0)THEN                                   ! if  "KGRID=1" selected
180
181 !
182 !*       2.1   Array allocation when called for the first time
183 !
184 !   1D Arrays
185 !
186   ALLOCATE(ZXTEM(1:IIU),STAT=OKZXTEM) ! RINT *,' OKZXTEM',OKZXTEM
187   ALLOCATE(ZYTEM(1:IJU),STAT=OKZYTEM) !PRINT *,' OKZYTEM',OKZYTEM
188   ALLOCATE(ZZTEM(1:IKU),STAT=OKZZTEM) !PRINT *,' OKZZTEM',OKZZTEM
189
190   ALLOCATE(ZDXTEM(1:IIU),STAT=OKZDXTEM) !PRINT *,' OKZDXTEM',OKZDXTEM
191   ALLOCATE(ZDYTEM(1:IJU),STAT=OKZDYTEM) !PRINT *,' OKZDYTEM',OKZDYTEM
192   ALLOCATE(ZDZTEM(1:IKU),STAT=OKZDZTEM) !PRINT *,' OKZDZTEM',OKZDZTEM
193
194   ALLOCATE(ZXMASS(1:IIU),STAT=OKZXMASS) !PRINT *,' OKZXMASS',OKZXMASS
195   ALLOCATE(ZYMASS(1:IJU),STAT=OKZYMASS) !PRINT *,' OKZYMASS',OKZYMASS
196   ALLOCATE(ZZMASS(1:IKU),STAT=OKZZMASS) !PRINT *,' OKZZMASS',OKZZMASS
197
198 !   2D Arrays
199 !
200   ALLOCATE(ZZSTEM(1:IIU,1:IJU),STAT=OKZZSTEM) !PRINT *,' OKZZSTEM',OKZZSTEM
201 !
202   IF(ALLOCATED(ZSCOEF))THEN
203     DEALLOCATE(ZSCOEF)
204   END IF
205     ALLOCATE(ZSCOEF(1:IIU,1:IJU),STAT=OKZSCOEF) !PRINT *,' OKZSCOEF',OKZSCOEF
206   IF(ALLOCATED(XXX))THEN
207     DEALLOCATE(XXX)
208   END IF
209     ALLOCATE(XXX(1:IIU,7))
210   IF(ALLOCATED(XXY))THEN
211     DEALLOCATE(XXY)
212   END IF
213     ALLOCATE(XXY(1:IJU,7))
214   IF(ALLOCATED(XXZ))THEN
215     DEALLOCATE(XXZ)
216   END IF
217     ALLOCATE(XXZ(1:IKU,7))
218   IF(ALLOCATED(XXDXHAT))THEN
219     DEALLOCATE(XXDXHAT)
220   END IF
221     ALLOCATE(XXDXHAT(1:IIU,7))
222   IF(ALLOCATED(XXDYHAT))THEN
223     DEALLOCATE(XXDYHAT)
224   END IF
225     ALLOCATE(XXDYHAT(1:IJU,7))
226
227 !   3D Arrays
228 !
229   IF(ALLOCATED(XXZS))THEN
230     DEALLOCATE(XXZS)
231   END IF
232     ALLOCATE(XXZS(1:IIU,1:IJU,7),STAT=OKXXZS) !PRINT *,' OKXXZS',OKXXZS
233   IF(ALLOCATED(ZSZ))THEN
234     DEALLOCATE(ZSZ)
235   END IF
236     ALLOCATE(ZSZ(1:IIU,1:IJU,IKU),STAT=OKZSZ) !PRINT *,' OKZSZ',OKZSZ
237 !
238 !*       2.2   Computes true altitudes on the W grid (KGRID=4)
239 !
240
241 IF(CSTORAGE_TYPE /= 'PG' .AND. CSTORAGE_TYPE /='SU')THEN
242   print *,' ******* COMPCOORD_FORDIACHRO  ZHAT(IKE+1) ',XZHAT(IKE+1)
243   CALL VERT_COORD(LSLEVE,XZS,XZSMT,XLEN1,XLEN2,XZHAT,ZSZ)
244 ENDIF
245 !
246 !*       2.3    Interpolates XHAT, YHAT, ZHAT at mass gridpoints
247 !
248
249   ZXMASS(1:IIU-1)=.5*(XXHAT(2:IIU)+XXHAT(1:IIU-1))
250   ZXMASS(IIU)=2.*ZXMASS(IIU-1)-ZXMASS(IIU-2)
251   ZYMASS(1:IJU-1)=.5*(XYHAT(2:IJU)+XYHAT(1:IJU-1))
252   ZYMASS(IJU)=2.*ZYMASS(IJU-1)-ZYMASS(IJU-2)
253   IF(IKU == 1)THEN
254     ZZMASS(1:IKU)=XZHAT(1:IKU)
255     !ZZMASS(1:IKU)=.5*(XZHAT(2:IKU)+XZHAT(1:IKU-1))  !!! size(XZHAT)=1 !!!
256   ELSE
257     ZZMASS(1:IKU-1)=.5*(XZHAT(2:IKU)+XZHAT(1:IKU-1))
258     ZZMASS(IKU)=2.*ZZMASS(IKU-1)-ZZMASS(IKU-2)
259   ENDIF
260
261 !
262 !*       2.4    Interpolates X, Y, Z, meshsizes, and topography 
263 !*              for all the KGRID selection  locations
264 !
265   DO JGRIDLOOP=1,7
266
267     SELECT CASE(JGRIDLOOP)
268    
269       CASE(1)
270         ZXTEM(:)=ZXMASS(:)
271         ZYTEM(:)=ZYMASS(:)
272         ZZTEM(:)=ZZMASS(:)
273         ZZSTEM(:,:)=XZS(:,:)
274   
275       CASE(2)
276         ZXTEM(:)=XXHAT(:)
277         ZYTEM(:)=ZYMASS(:)
278         ZZTEM(:)=ZZMASS(:)
279         ZZSTEM(2:IIU,:)=.5*(XZS(2:IIU,:)+XZS(1:IIU-1,:))
280         ZZSTEM(1,:)=XZS(1,:)
281   
282       CASE(3)
283         ZXTEM(:)=ZXMASS(:)
284         ZYTEM(:)=XYHAT(:)
285         ZZTEM(:)=ZZMASS(:)
286         ZZSTEM(:,2:IJU)=.5*(XZS(:,2:IJU)+XZS(:,1:IJU-1))
287         ZZSTEM(:,1)=XZS(:,1)
288    
289       CASE(4)
290         ZXTEM(:)=ZXMASS(:)
291         ZYTEM(:)=ZYMASS(:)
292         ZZTEM(:)=XZHAT(:)
293         ZZSTEM(:,:)=XZS(:,:)
294    
295       CASE(5)
296         ZXTEM(:)=XXHAT(:)
297         ZYTEM(:)=XYHAT(:)
298         ZZTEM(:)=ZZMASS(:)
299         ZZSTEM(2:IIU,:)=.5*(XZS(2:IIU,:)+XZS(1:IIU-1,:))
300         ZZSTEM(1,:)=XZS(1,:)
301         ZZSTEM(:,2:IJU)=.5*(ZZSTEM(:,2:IJU)+ZZSTEM(:,1:IJU-1))
302         ZZSTEM(:,1)=ZZSTEM(:,2)
303   
304       CASE(6)
305         ZXTEM(:)=XXHAT(:)
306         ZYTEM(:)=ZYMASS(:)
307         ZZTEM(:)=XZHAT(:)
308         ZZSTEM(2:IIU,:)=.5*(XZS(2:IIU,:)+XZS(1:IIU-1,:))
309         ZZSTEM(1,:)=XZS(1,:)
310   
311       CASE(7)
312         ZXTEM(:)=ZXMASS(:)
313         ZYTEM(:)=XYHAT(:)
314         ZZTEM(:)=XZHAT(:)
315         ZZSTEM(:,2:IJU)=.5*(XZS(:,2:IJU)+XZS(:,1:IJU-1))
316         ZZSTEM(:,1)=XZS(:,1)
317   
318     END SELECT
319    
320     ZDXTEM(1:IIU-1)=ZXTEM(2:IIU)-ZXTEM(1:IIU-1)
321 !
322 ! NOTICE: An extra meshlength is added to the max size of the arrays
323 ! in order to avoid a lot of testing hereafter...
324 !
325     ZDXTEM(IIU)=ZDXTEM(IIU-1)
326     ZDYTEM(1:IJU-1)=ZYTEM(2:IJU)-ZYTEM(1:IJU-1)
327     ZDYTEM(IJU)=ZDYTEM(IJU-1)
328     IF(IKU /= 1)THEN
329     ZDZTEM(1:IKU-1)=ZZTEM(2:IKU)-ZZTEM(1:IKU-1)
330     ZDZTEM(IKU)=ZDZTEM(IKU-1)
331     ENDIF
332   
333 ! X, Y, Z as functions of KGRID
334     XXX(:,JGRIDLOOP)=ZXTEM
335     XXY(:,JGRIDLOOP)=ZYTEM
336     XXZ(:,JGRIDLOOP)=ZZTEM
337   
338 ! Topography as a function of KGRID
339     XXZS(:,:,JGRIDLOOP)=ZZSTEM
340    
341 ! Meshsizes as functions of KGRID
342     XXDXHAT(:,JGRIDLOOP)=ZDXTEM(:)
343     XXDYHAT(:,JGRIDLOOP)=ZDYTEM(:)
344   
345   ENDDO
346    
347   DEALLOCATE(ZXMASS,ZYMASS,ZZMASS)
348   DEALLOCATE(ZXTEM,ZYTEM,ZZTEM)
349   DEALLOCATE(ZDXTEM,ZDYTEM,ZDZTEM)
350   DEALLOCATE(ZZSTEM)
351
352 !-------------------------------------------------------------------------------
353 !
354 !*       3.   CALCULATIONS PERFORMED FOR ALL SUBSEQUENT CALLS
355 !             -----------------------------------------------
356 !
357 ELSE                                            ! else if KGRID =/=1 selected
358    
359 ! True altitudes
360
361   SELECT CASE(KGRID)
362    
363     CASE(1)
364       XZZ(:,:,1:IKU-1)=0.5*(ZSZ(:,:,1:IKU-1)+ZSZ(:,:,2:IKU))
365       XZZ(:,:,IKU)=2.*XZZ(:,:,IKU-1)-XZZ(:,:,IKU-2)
366    
367     CASE(2)
368       XZZ(:,:,1:IKU-1)=0.5*(ZSZ(:,:,1:IKU-1)+ZSZ(:,:,2:IKU))
369       XZZ(:,:,IKU)=2.*XZZ(:,:,IKU-1)-XZZ(:,:,IKU-2)
370       XZZ(2:IIU,:,:)=0.5*(XZZ(2:IIU,:,:)+XZZ(1:IIU-1,:,:))
371       XZZ(1,:,:)=2*XZZ(2,:,:)-XZZ(3,:,:)
372
373     CASE(3)
374       XZZ(:,:,1:IKU-1)=0.5*(ZSZ(:,:,1:IKU-1)+ZSZ(:,:,2:IKU))
375       XZZ(:,:,IKU)=2.*XZZ(:,:,IKU-1)-XZZ(:,:,IKU-2)
376       XZZ(:,2:IJU,:)=0.5*(XZZ(:,2:IJU,:)+XZZ(:,1:IJU-1,:))
377       XZZ(:,1,:)=2*XZZ(:,2,:)-XZZ(:,3,:)
378
379     CASE(4)
380       XZZ(:,:,:)=ZSZ(:,:,:)
381
382     CASE(5)
383       XZZ(:,:,1:IKU-1)=0.5*(ZSZ(:,:,1:IKU-1)+ZSZ(:,:,2:IKU))
384       XZZ(:,:,IKU)=2.*XZZ(:,:,IKU-1)-XZZ(:,:,IKU-2)
385       XZZ(2:IIU,:,:)=0.5*(XZZ(2:IIU,:,:)+XZZ(1:IIU-1,:,:))
386       XZZ(1,:,:)=2*XZZ(2,:,:)-XZZ(3,:,:)
387       XZZ(:,2:IJU,:)=0.5*(XZZ(:,2:IJU,:)+XZZ(:,1:IJU-1,:))
388       XZZ(:,1,:)=2*XZZ(:,2,:)-XZZ(:,3,:)
389
390     CASE(6)
391       XZZ(2:IIU,:,:)=0.5*(ZSZ(2:IIU,:,:)+ZSZ(1:IIU-1,:,:))
392       XZZ(1,:,:)=2*XZZ(2,:,:)-XZZ(3,:,:)
393
394     CASE(7)
395       XZZ(:,2:IJU,:)=0.5*(ZSZ(:,2:IJU,:)+ZSZ(:,1:IJU-1,:))
396       XZZ(:,1,:)=2*XZZ(:,2,:)-XZZ(:,3,:)
397
398   END SELECT
399
400 END IF                                                ! End KGRID selection
401 !
402 !---------------------------------------------------------------------------
403 !
404 !*       4.   EXIT
405 !             ----
406 !
407 RETURN
408 END SUBROUTINE COMPCOORD_FORDIACHRO