G.Delautier 19/09/2017 : bug MNH2LPDM
[MNH-git_open_source-lfs.git] / src / MNH / mnh2lpdm_ini.f90
1 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
3 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
4 !MNH_LIC for details. version 1.
5 !     ######spl
6         SUBROUTINE MNH2LPDM_INI(HFM1,HFM2,HFLOG,KFLOG,KFGRI,KFDAT,KVERB)
7 !--------------------------------------------------------------------------
8 !* MNH2S2_INI    : INITIALISATION DU COUPLAGE MESO-NH / LPDM.
9 !
10 !  Auteur        : Francois BONNARDOT, DP/SERV/ENV
11 !  Creation      : 04.01.2009 (mnh2s2_ini.f90)
12 !
13 !
14 !       Arguments explicites.
15 !       ---------------------
16 !       HFM1,HFM2     Nom du premier et du dernier Fichier FM à traiter.
17 !       HFLOG   Nom du fichier-log.
18 !       KFLOG   Unite logique du fichier-log.
19 !       IFGRI   Unite logique du fichier GRILLE.
20 !       IFDAT   Unite logique du fichier DATE
21 !       KVERB   Niveau de blabla.
22 !
23 !--------------------------------------------------------------------------
24 !
25 !
26 !
27 !*      0.  INITIALISATION.
28 !           ---------------
29 !
30 !*      0.1 Modules.
31 !
32 USE MODE_MODELN_HANDLER
33 USE MODD_MNH2LPDM
34 !
35 USE MODE_FM
36 USE MODE_IO_ll
37 !
38 USE MODD_CST
39 USE MODD_TIME
40 USE MODD_PARAMETERS
41 USE MODD_GRID
42 USE MODD_DIM_n
43 USE MODD_GRID_n
44 USE MODD_TIME_n
45 USE MODD_LUNIT
46 !
47 USE MODI_INI_CST
48 USE MODI_READ_HGRID
49 USE MODE_FMREAD
50 !
51 USE MODE_GRIDPROJ
52 USE MODI_XYTOLATLON
53 !
54 USE MODI_TEMPORAL_DIST
55 !
56 !*      0.2 Arguments.
57 !
58 IMPLICIT NONE
59 !
60 CHARACTER(LEN=*), INTENT(IN)  :: HFM1,HFM2,HFLOG
61 INTEGER,          INTENT(IN)  :: KFLOG,KFGRI,KFDAT,KVERB
62 !
63 !
64 !*      0.3 Variables locales.
65 !
66 CHARACTER(LEN=28)     :: YNAME,YDAD     ! Noms du FM et de son papa.
67 CHARACTER(LEN=100)    :: YCOM           ! Commentaire MNH.
68 CHARACTER(LEN=2)      :: YSTORAGE       ! Type de variable.
69 !
70 INTEGER, DIMENSION(3) :: IDATMDL,IDATCUR1,IDATCUR2 ! Date  du modele, Date courante
71 REAL                  :: ZSECMDL,ZSECCUR1,ZSECCUR2 ! Heure du modele, Heure courante
72 REAL                  :: ZECHEANCE1,ZECHEANCE2     ! dist temp date modele - date courante
73 INTEGER               :: IHHMDL,IMNMDL,ISSMDL ! h - mn - s du model
74 INTEGER               :: IHHCUR1,IMNCUR1,ISSCUR1
75 INTEGER               :: IHHCUR2,IMNCUR2,ISSCUR2
76 CHARACTER(LEN=14)     :: YDATMDL,YDATCUR1,YDATCUR2
77 !
78 INTEGER               :: INBART,IGRID,ILONCOM,IREP
79 REAL                  :: XLATOR,XLONOR,XPTLAT,XPTLON
80 REAL                  :: XXPTSOMNH,XYPTSOMNH
81 INTEGER               :: JI,JJ,JK,a
82 INTEGER               :: b,c,I
83 INTEGER, DIMENSION(:,:), ALLOCATABLE   :: TAB2D
84 INTEGER, DIMENSION(:), ALLOCATABLE   :: TAB1D
85 !
86 !
87 !
88 !*      1.  INITIALISATION.
89 !           ---------------
90 !
91 CALL INI_CST
92 !
93 CALL GOTO_MODEL(1)
94 !
95 !
96 !*      2.  DONNEES MESO-NH.
97 !           ----------------
98 !
99 !*      2.1 Ouverture du fichier Meso-NH.
100 !
101 CALL FMOPEN_LL(HFM1,'READ',HFLOG,0,2,KVERB,INBART,IREP)
102 CALL FMOPEN_LL(HFM2,'READ',HFLOG,0,2,KVERB,INBART,IREP)
103 !
104 !
105 !*      2.2 Date et heure du modele.
106 !
107 CALL FMREAD(HFM1,'DTEXP%TDATE',HFLOG,'--',IDATMDL,IGRID,ILONCOM,YCOM,IREP)
108 CALL FMREAD(HFM1,'DTEXP%TIME', HFLOG,'--',ZSECMDL,IGRID,ILONCOM,YCOM,IREP)
109 CALL FMREAD(HFM1,'DTCUR%TDATE',HFLOG,'--',IDATCUR1,IGRID,ILONCOM,YCOM,IREP)
110 CALL FMREAD(HFM1,'DTCUR%TIME', HFLOG,'--',ZSECCUR1,IGRID,ILONCOM,YCOM,IREP)
111 CALL FMREAD(HFM2,'DTCUR%TDATE',HFLOG,'--',IDATCUR2,IGRID,ILONCOM,YCOM,IREP)
112 CALL FMREAD(HFM2,'DTCUR%TIME', HFLOG,'--',ZSECCUR2,IGRID,ILONCOM,YCOM,IREP)
113 !
114 CALL TEMPORAL_DIST(IDATCUR1(1),IDATCUR1(2),IDATCUR1(3), &
115                    ZSECCUR1,IDATMDL(1),IDATMDL(2),IDATMDL(3), &
116                    ZSECMDL,ZECHEANCE1)
117 CALL TEMPORAL_DIST(IDATCUR2(1),IDATCUR2(2),IDATCUR2(3), &
118                    ZSECCUR2,IDATMDL(1),IDATMDL(2),IDATMDL(3), &
119                    ZSECMDL,ZECHEANCE2)
120 !
121 IHHMDL=INT(ZSECMDL/3600)
122 IMNMDL=INT((ZSECMDL-IHHMDL*3600)/60)
123 ISSMDL=INT(ZSECMDL-IHHMDL*3600-IMNMDL*60)
124 IHHCUR1=INT(ZSECCUR1/3600)
125 IMNCUR1=INT((ZSECCUR1-IHHCUR1*3600)/60)
126 ISSCUR1=INT(ZSECCUR1-IHHCUR1*3600-IMNCUR1*60)
127 IHHCUR2=INT(ZSECCUR2/3600)
128 IMNCUR2=INT((ZSECCUR2-IHHCUR2*3600)/60)
129 ISSCUR2=INT(ZSECCUR2-IHHCUR2*3600-IMNCUR2*60)
130 !
131 WRITE(YDATMDL,'(I4.4,5I2.2)') IDATMDL(1), IDATMDL(2), IDATMDL(3), &
132                                IHHMDL, IMNMDL, ISSMDL
133 WRITE(YDATCUR1,'(I4.4,5I2.2)') IDATCUR1(1), IDATCUR1(2), IDATCUR1(3), &
134                                IHHCUR1, IMNCUR1, ISSCUR1
135 WRITE(YDATCUR2,'(I4.4,5I2.2)') IDATCUR2(1), IDATCUR2(2), IDATCUR2(3), &
136                                IHHCUR2, IMNCUR2, ISSCUR2
137
138 NMDLAA=MOD( IDATMDL(1), 100 )  ! Annee arrondi a 2 chiffres.
139 NMDLMM=IDATMDL(2)
140 NMDLJJ=IDATMDL(3)
141 NMDLSS=NINT(ZSECMDL)
142 !
143 !*      Heure du modele arrondie a 5 minutes pres.
144 !
145 NMDLMN = NINT( (FLOAT(NMDLSS)/60.0)/5.0 )*5
146 NMDLSS = 0
147 NMDLHH =NMDLMN/60
148 NMDLMN =NMDLMN-NMDLHH*60
149 !
150 !*      2.3 Grille horizontale.
151 !
152 CALL READ_HGRID(1,HFM1,YNAME,YDAD,YSTORAGE)
153 IF (YNAME == YDAD) THEN
154 IGRILLE=1
155 ELSE
156 IGRILLE=2
157 ENDIF
158 print*,IGRILLE
159 !
160 ! Lecture grille horizontale
161 !
162 NIU=NIMAX+2*JPHEXT
163 NJU=NJMAX+2*JPHEXT
164 NIB=1+JPHEXT
165 NJB=1+JPHEXT
166 NIE=NIU-JPHEXT
167 NJE=NJU-JPHEXT
168 !
169 !
170 !*      2.4 Nombre de niveaux-verticaux.
171 !
172 CALL FMREAD(HFM1,'KMAX',HFLOG,'--',NKMAX,IGRID,ILONCOM,YCOM,IREP)
173 !WRITE(KFLOG,*) '%%% MNH2S2_INI Lecture du nombre de niveau OK.'
174 !
175 NKU = NKMAX+2*JPVEXT
176 NKB = 1+JPVEXT
177 NKE = NKU-JPVEXT
178 !
179 !
180 !*      2.5 Allocations Meso-NH.
181 !
182 ALLOCATE( XZHAT(NKU) )
183 ALLOCATE( XZS(NIU,NJU) )
184 ALLOCATE( XZ0(NIU,NJU) )
185 ALLOCATE( XUT(NIU,NJU,NKU))
186 ALLOCATE( XVT(NIU,NJU,NKU))
187 ALLOCATE( XWT(NIU,NJU,NKU))
188 ALLOCATE( XTHT(NIU,NJU,NKU))
189 ALLOCATE( XTKET(NIU,NJU,NKU))
190 ALLOCATE( XLM(NIU,NJU,NKU))
191 ALLOCATE( XDISSIP(NIU,NJU,NKU))
192 ALLOCATE( XWPTHP(NIU,NJU,NKU))
193 ALLOCATE( XRMVT(NIU,NJU,NKU))
194 ALLOCATE( XRMCT(NIU,NJU,NKU))
195 ALLOCATE( XRMRT(NIU,NJU,NKU))
196 ALLOCATE( XINRT(NIU,NJU))
197 ALLOCATE( XSFU(NIU,NJU))
198 ALLOCATE( XSFV(NIU,NJU))
199 !
200 !*      2.6 Decoupage vertical.
201 !
202 CALL FMREAD(HFM1,'ZHAT',HFLOG,'--',XZHAT,IGRID,ILONCOM,YCOM,IREP)
203 !
204 !*      2.7 Orographie. 
205 !
206 CALL FMREAD(HFM1,'ZS',HFLOG,'XY',XZS,IGRID,ILONCOM,YCOM,IREP)
207 !
208 !*      2.8 Rugosite Z0. 
209 !
210 CALL FMREAD(HFM1,'Z0',HFLOG,'XY',XZ0,IGRID,ILONCOM,YCOM,IREP)
211 !
212 XXPTSOMNH=XXHAT(1)+(XXHAT(2)-XXHAT(1))/2
213 XYPTSOMNH=XYHAT(1)+(XYHAT(2)-XYHAT(1))/2
214 CALL SM_LATLON(XLATORI,XLONORI,XXPTSOMNH,XYPTSOMNH,XLATOR,XLONOR)
215 !
216 !*      2.9  DOMAINE D'EXTRACTION.
217 !           ---------------------
218 !
219 NSIB   = NIB
220 NSIE   = NIE
221 NSJB   = NJB
222 NSJE   = NJE
223 !
224 NSIMAX = NSIE-NSIB+1
225 NSJMAX = NSJE-NSJB+1
226 !
227 !
228 !*      3. Impression de controle Meso-NH.
229 !           -------------------------------
230 !
231 !           Domaine horizontal Meso-NH.
232 !modif 12.2014 : passage a 1 seul domaine MesoNH
233 !           ---------------------------
234 WRITE(KFLOG,'(I1,a12)') IGRILLE,'      ngrid '
235 !WRITE(KFLOG,'(a13)') '2      ngrids'
236 WRITE(KFLOG,'(a13)') '1      ngrids'
237 WRITE(KFLOG,'(i4,3x,a6)') NSIMAX,'nx    '
238 WRITE(KFLOG,'(i4,3x,a6)') NSJMAX,'ny    '
239 WRITE(KFLOG,'(i4,3x,a6)') NKU-2,'nz    '
240 WRITE(KFLOG,'(i4,3x,a6)') NKU-3,'nzg   ' 
241 WRITE(KFLOG,'(a13)') '12     npatch'
242 WRITE(KFLOG,'(a13)') '0      icloud'
243 WRITE(KFLOG,'(a11)') '0.0  wlon  '
244 WRITE(KFLOG,'(a11)') '45.0 rnlat '
245 WRITE(KFLOG,'(f10.1,3x,a6)') XZHAT(NKE),'s     '
246 WRITE(KFLOG,'(f8.0,a8)') ZECHEANCE1,'  time1 '
247 WRITE(KFLOG,'(f8.0,a8)') ZECHEANCE2,'  time2 '
248 WRITE(KFLOG,'(a13)') '3600    dtmet '
249 WRITE(KFLOG,'(a13)') 'm       tunits'
250 WRITE(KFLOG,'(a13)') '12     nvout '
251 WRITE(KFLOG,'(6x,a8,i4)') 'u       ',1
252 WRITE(KFLOG,'(6x,a8,i4)') 'v       ',1
253 WRITE(KFLOG,'(6x,a8,i4)') 'w       ',1
254 WRITE(KFLOG,'(6x,a8,i4)') 'tp      ',1
255 WRITE(KFLOG,'(6x,a8,i4)') 'tke     ',1
256 WRITE(KFLOG,'(6x,a8,i4)') 'uu      ',1
257 WRITE(KFLOG,'(6x,a8,i4)') 'vv      ',1
258 WRITE(KFLOG,'(6x,a8,i4)') 'ww      ',1
259 WRITE(KFLOG,'(6x,a8,i4)') 'tlx     ',1
260 WRITE(KFLOG,'(6x,a8,i4)') 'tly     ',1
261 WRITE(KFLOG,'(6x,a8,i4)') 'tlz     ',1
262 WRITE(KFLOG,'(6x,a8,i4)') 'intopr  ',1
263 WRITE(KFLOG,*) '  grid structure'
264 !
265 !*      4.  FICHIER METEO.
266 !           --------------
267 !
268 !*      4.1 Allocations.
269 !
270 ALLOCATE( XSHAUT(NKMAX))
271 ALLOCATE( XSREL(NSIMAX,NSJMAX) )
272 ALLOCATE( XSZ0(NSIMAX,NSJMAX) )
273 ALLOCATE( XSCORIOZ (NSIMAX,NSJMAX) )
274 ALLOCATE( XSPHI(NSIMAX,NSJMAX,NKMAX) )
275 ALLOCATE( XSU(NSIMAX,NSJMAX,NKMAX))
276 ALLOCATE( XSV(NSIMAX,NSJMAX,NKMAX))
277 ALLOCATE( XSW(NSIMAX,NSJMAX,NKMAX))
278 ALLOCATE( XSTH(NSIMAX,NSJMAX,NKMAX))
279 ALLOCATE( XSTKE(NSIMAX,NSJMAX,NKMAX))
280 ALLOCATE( XSLM(NSIMAX,NSJMAX,NKMAX))
281 ALLOCATE( XSDISSIP(NSIMAX,NSJMAX,NKMAX))
282 ALLOCATE( XSWPTHP(NSIMAX,NSJMAX,NKMAX))
283 ALLOCATE( XSRMV(NSIMAX,NSJMAX,NKMAX))
284 ALLOCATE( XSRMC(NSIMAX,NSJMAX,NKMAX))
285 ALLOCATE( XSRMR(NSIMAX,NSJMAX,NKMAX))
286 ALLOCATE( XSINRT(NSIMAX,NSJMAX))
287 ALLOCATE( XSSFU(NSIMAX,NSJMAX))
288 ALLOCATE( XSSFV(NSIMAX,NSJMAX))
289 ALLOCATE( XSTIMEW(NSIMAX,NSJMAX,NKMAX))
290 ALLOCATE( XSTIMEU(NSIMAX,NSJMAX,NKMAX))
291 ALLOCATE( XSSIGW(NSIMAX,NSJMAX,NKMAX))
292 ALLOCATE( XSSIGU(NSIMAX,NSJMAX,NKMAX))
293 ALLOCATE( XSUSTAR(NSIMAX,NSJMAX))
294 ALLOCATE( XSWSTAR(NSIMAX,NSJMAX))
295 ALLOCATE( XSHMIX(NSIMAX,NSJMAX))
296 ALLOCATE( XSLMO(NSIMAX,NSJMAX))
297 ALLOCATE( XSTHETAV(NKMAX))
298
299 !
300 !   4.2.    Nombre de niveaux en Z
301 !
302 XSHAUT(1:NKMAX) = (XZHAT(NKB:NKE)+XZHAT(NKB+1:NKE+1))/2.
303 print*,"niveaux hauteur"
304 DO JK=1,NKMAX
305 print*,XSHAUT(JK)
306 ENDDO
307 !
308 !   4.3.    Calcul du tableau contenant les coef. de coriolis de la grille
309 !
310 DO JI=NSIB,NSIE ; DO JJ=NSJB,NSJE
311    CALL SM_LATLON(XLATORI,XLONORI,XXHAT(JI),XYHAT(JJ),XPTLAT,XPTLON)
312    XSCORIOZ(JI-1,JJ-1)=2.*XOMEGA*SIN(XPTLAT*XPI/180.)
313 ENDDO ; ENDDO
314 !
315 !
316 !* 4.4 Geometrie de la grille et positionnement.
317 !
318 !
319 ! On a besoin du point sud-ouest, c'est-a-dire de l'angle inferieur gauche
320 ! du domaine physique de la maille "en bas a gauche". Ca tombe bien, on
321 ! va travailler avec les XXHAT et les XYHAT directement.
322 !
323 XPASXM = XXHAT(2)-XXHAT(1)     ! Pas selon X en metres.
324 XPASYM = XYHAT(2)-XYHAT(1)     ! Pas selon Y en metres.
325 ZMAILLE = MAX(XPASXM,XPASYM)
326 !
327 !* 4.5 Constantes et champs constants.
328 !
329 !* Relief.
330 !
331 XSREL(:,:)  =  XZS(NSIB:NSIE,NSJB:NSJE)
332 !
333 !* Geopotentiel PHI
334 !
335 print*,"Geopotentiel"
336 DO JK=1,NKMAX
337 XSPHI(:,:,JK) = (XSREL(:,:)+XSHAUT(JK))*XG
338 print*,MINVAL(XSPHI(:,:,JK)),MAXVAL(XSPHI(:,:,JK))
339 ENDDO
340 !
341 !* Rugosite.
342 !
343 XSZ0(:,:)  =  XZ0(NSIB:NSIE,NSJB:NSJE)
344 print*,"Rugosite"
345 print*,MINVAL(XSZ0),MAXVAL(XSZ0)
346 !
347 !* 5   FICHIER DATES.
348 !      -------------
349 !
350 WRITE(KFDAT,'(A14)') YDATMDL 
351 WRITE(KFDAT,'(A14)') YDATCUR1 
352 WRITE(KFDAT,'(A14)') YDATCUR2
353 !
354 !* 5.  FICHIER GRILLE.
355 !      --------------
356 !
357 !
358 !* 5.1 Infos franchement utiles.
359 !
360 WRITE(KFGRI,'(F15.8,1X,A)') &
361                   XLON0,   'XLON0   Longitude reference (deg.deci.)'
362 WRITE(KFGRI,'(F15.8,1X,A)') &
363                   XLAT0,   'XLAT0   Latitude  reference (deg.deci.)'
364 WRITE(KFGRI,'(F15.8,1X,A)') &
365                   XBETA,   'XBETA   Rotation  grille    (deg.deci.)'
366 WRITE(KFGRI,'(F15.8,1X,A)') XRPK,    'XRPK    Facteur de conicite'
367 WRITE(KFGRI,'(F15.8,1X,A)') &
368                   XLONOR,  'XLONOR  Longitude origine   (deg.deci.)'
369 WRITE(KFGRI,'(F15.8,1X,A)') &
370                   XLATOR,  'XLATOR  Latitude  origine   (deg.deci.)'
371 WRITE(KFGRI,'(F15.1,1X,A)') XXHAT(1),'XHAT(1) Coord. Cartesienne  (m)'
372 WRITE(KFGRI,'(F15.1,1X,A)') XXHAT(2),'XHAT(2) Coord. Cartesienne  (m)'
373 WRITE(KFGRI,'(F15.1,1X,A)') XYHAT(1),'YHAT(1) Coord. Cartesienne  (m)'
374 WRITE(KFGRI,'(F15.1,1X,A)') XYHAT(2),'YHAT(2) Coord. Cartesienne  (m)'
375 !
376 print*,"GRILLE"
377 print*,"LON0 : ",XLON0
378 print*,"LAT0 : ",XLAT0
379 print*,"BETA : ",XBETA
380 print*,"RPK  : ",XRPK
381 print*,"LONOR: ",XLONOR
382 print*,"LATOR: ",XLATOR
383 !
384 !* 5.2  Points de grille x y z zg
385 !
386 WRITE(KFLOG,*)NSIMAX,' gridpoints in x direction'
387 WRITE(KFLOG,'(8f10.0)')XXHAT(NSIB:NSIE)
388 WRITE(KFLOG,*)NSJMAX,' gridpoints y direction'
389 WRITE(KFLOG,'(8f10.0)')XYHAT(NSJB:NSJE)
390 WRITE(KFLOG,*)NKMAX,'  main gridpoints in z direction'
391 WRITE(KFLOG,'(8f10.2)')XSHAUT(1:NKMAX)
392 WRITE(KFLOG,'(i4,3x,a38)')NKU-2,'intermediate gridpoints in z direction'
393 WRITE(KFLOG,'(8f10.2)')XZHAT(2:NKU-1)
394 WRITE(KFLOG,*)'   =================================================='
395 !
396 !     Topographie
397 !
398 WRITE(KFLOG,*) 'TERRAIN TOPOGRAPHY'
399 c=1
400 a=0
401 !modif 12/2014 : passage a une grille haute resolution MesoNH, on depasse 99
402 !300 format(i2,'|',18i4)
403 300 format(i3,'|',18i5)
404 !400 format(i2,'|',18(f4.2))
405 !400 format(i3,'|',18(f5.2))
406 !301 format(3x,18('__',i2))
407 301 format(3x,18('__',i3))
408 ALLOCATE(TAB2D(NSIMAX,NSJMAX))
409 ALLOCATE(TAB1D(NSIMAX))
410 DO I=1,NSIMAX
411    TAB1D(I)=I
412 ENDDO
413 TAB2D(:,:) = NINT(XSREL(:,:))
414 DO WHILE (c.lt.(NSIMAX+1))
415   DO b=NSJB,NSJE
416      IF ((c+17).LT.(NSIMAX+1)) then
417         a=NSJMAX-b+NSJB
418         WRITE(KFLOG,300) a,TAB2D(c:c+17,a)
419      ELSE  
420         a=NSJMAX-b+NSJB
421         WRITE(KFLOG,300) a,TAB2D(c:NSIMAX,a)
422      ENDIF
423   ENDDO
424 IF ((c+17).LT.(NSIMAX+1)) then
425    WRITE(KFLOG,301) TAB1D(c:c+17)
426 ELSE
427    WRITE(KFLOG,301) TAB1D(c:NSIMAX)
428 ENDIF
429
430 c=c+18
431 ENDDO
432 !
433 DEALLOCATE(TAB2D)
434 DEALLOCATE(TAB1D)
435 DEALLOCATE(XZS)
436 DEALLOCATE(XZ0)
437 DEALLOCATE(XZHAT)
438 !
439 !        Fermeture du fichier Meso-NH.
440 !
441 CALL FMCLOS_LL(HFM1,'KEEP',HFLOG,IREP)
442 CALL FMCLOS_LL(HFM2,'KEEP',HFLOG,IREP)
443 !
444 !
445 !-------------------------------------------'
446 print*,'          FIN MNH2LPDM_INI'
447 !-------------------------------------------'
448 !
449 !
450 END SUBROUTINE MNH2LPDM_INI