SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
mode_gridtype_ign.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 ! ##############################
6  MODULE mode_gridtype_ign
7 ! ##############################
8 !
9 !############################################################################
10 !############################################################################
11 !############################################################################
12 !
13 USE yomhook ,ONLY : lhook, dr_hook
14 USE parkind1 ,ONLY : jprb
15 !
16  CONTAINS
17 !############################################################################
18 !############################################################################
19 !############################################################################
20 ! ####################################################################
21  SUBROUTINE put_gridtype_ign(PGRID_PAR,KLAMBERT,PX,PY,PDX,PDY, &
22  kdimx,kdimy,pxall,pyall )
23 ! ####################################################################
24 !
25 !!**** *PUT_GRIDTYPE_IGN* - routine to store in PGRID_PAR the horizontal grid
26 !!
27 !! AUTHOR
28 !! ------
29 !! E. Martin *Meteo France*
30 !!
31 !! MODIFICATIONS
32 !! -------------
33 !! Original 10/2007
34 !! 02/2011 Correction de la longitude d'origine pourle cas L93 (A. Lemonsu)
35 !! 07/2011 add maximum domain dimension for output (B. Decharme)
36 ! 01/2016 Correction de la valeur de l'excentricite pour L93 (V. Masson)
37 ! 01/2016 Correction de la valeur du rayon terrestre pour Lamberts 1 a 4 (V. Masson)
38 ! 01/2016 Correction d'une boucle pour parallelisation (V. Masson)
39 !-------------------------------------------------------------------------------
40 !
41 !* 0. DECLARATIONS
42 ! ------------
43 !
44 IMPLICIT NONE
45 !
46 !
47 !* 0.1 Declarations of arguments
48 ! -------------------------
49 !
50 INTEGER, INTENT(IN) :: klambert ! Lambert type
51  ! 1 : Lambert I
52  ! 2 : Lambert II
53  ! 3 : Lambert III
54  ! 4 : Lambert IV
55  ! 5 : Extended Lambert II
56  ! 6 : Lambert 93
57 REAL, DIMENSION(:), INTENT(IN) :: px ! X coordinate of grid mesh center
58 REAL, DIMENSION(:), INTENT(IN) :: py ! Y coordinate of grid mesh center
59 REAL, DIMENSION(:), INTENT(IN) :: pdx ! X grid mesh size
60 REAL, DIMENSION(:), INTENT(IN) :: pdy ! Y grid mesh size
61 INTEGER, INTENT(IN) :: kdimx ! maximum domain length in X
62 INTEGER, INTENT(IN) :: kdimy ! maximum domain length in Y
63 REAL, DIMENSION(KDIMX), INTENT(IN) :: pxall! maximum domain X coordinate of grid mesh
64 REAL, DIMENSION(KDIMY), INTENT(IN) :: pyall! maximum domain Y coordinate of grid mesh
65 REAL, DIMENSION(:), POINTER :: pgrid_par! parameters defining this grid
66 !
67 !
68 !* 0.2 Declarations of local variables
69 ! -------------------------------
70 !
71 INTEGER :: il ! number of points
72 REAL(KIND=JPRB) :: zhook_handle
73 !-------------------------------------------------------------------------------
74 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:PUT_GRIDTYPE_IGN',0,zhook_handle)
75 il = SIZE(px)
76 ALLOCATE(pgrid_par(4*il+4+kdimx+kdimy))
77 pgrid_par(1) = float(klambert)
78 pgrid_par(2) = float(il)
79 pgrid_par(3:il+2) = px(:)
80 pgrid_par(il+3:2*il+2) = py(:)
81 pgrid_par(2*il+3:3*il+2) = pdx(:)
82 pgrid_par(3*il+3:4*il+2) = pdy(:)
83 pgrid_par(4*il+3) = float(kdimx)
84 pgrid_par(4*il+4) = float(kdimy)
85 pgrid_par(4*il+5:4*il+4+kdimx) = pxall(:)
86 pgrid_par(4*il+5+kdimx:4*il+4+kdimx+kdimy) = pyall(:)
87 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:PUT_GRIDTYPE_IGN',1,zhook_handle)
88 !
89 !-------------------------------------------------------------------------------
90 END SUBROUTINE put_gridtype_ign
91 !############################################################################
92 !############################################################################
93 !############################################################################
94 ! ####################################################################
95  SUBROUTINE get_gridtype_ign(PGRID_PAR,KLAMBERT,KL,PX,PY,PDX,PDY,&
96  kdimx,kdimy,pxall,pyall )
97 ! ####################################################################
98 !
99 !!**** *GET_GRIDTYPE_IGN* - routine to get from PGRID_PAR the horizontal grid
100 !!
101 !! AUTHOR
102 !! ------
103 !! E. Martin *Meteo France*
104 !!
105 !! MODIFICATIONS
106 !! -------------
107 !! Original 10/2007
108 !! 07/2011 add maximum domain dimension for output (B. Decharme)
109 !-------------------------------------------------------------------------------
110 !
111 !* 0. DECLARATIONS
112 ! ------------
113 !
114 IMPLICIT NONE
115 !
116 !
117 !* 0.1 Declarations of arguments
118 ! -------------------------
119 !
120 INTEGER, INTENT(OUT), OPTIONAL :: kl ! number of points
121 INTEGER, INTENT(OUT), OPTIONAL :: klambert ! Lambert type
122  ! 1 : Lambert I
123  ! 2 : Lambert II
124  ! 3 : Lambert III
125  ! 4 : Lambert IV
126  ! 5 : Extended Lambert II
127  ! 6 : Lambert 93
128 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: px ! X coordinate of grid mesh center
129 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: py ! Y coordinate of grid mesh center
130 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: pdx ! X grid mesh size
131 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: pdy ! Y grid mesh size
132 INTEGER, INTENT(OUT), OPTIONAL :: kdimx ! maximum domain length in X
133 INTEGER, INTENT(OUT), OPTIONAL :: kdimy ! maximum domain length in Y
134 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: pxall ! maximum domain X coordinate of grid mesh
135 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: pyall ! maximum domain Y coordinate of grid mesh
136 
137 REAL, DIMENSION(:), INTENT(IN) :: pgrid_par! parameters defining this grid
138 !
139 !
140 !* 0.2 Declarations of local variables
141 ! -------------------------------
142 !
143 INTEGER :: il, idimx, idimy
144 !
145 REAL(KIND=JPRB) :: zhook_handle
146 !-------------------------------------------------------------------------------
147 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:GET_GRIDTYPE_IGN',0,zhook_handle)
148 IF (present(klambert)) klambert = nint(pgrid_par(1))
149 IF (present(kl)) kl = nint(pgrid_par(2))
150 !
151 il = nint(pgrid_par(2))
152 !
153 IF (present(px)) px(:) = pgrid_par(2+1:2+il)
154 !
155 IF (present(py)) py(:) = pgrid_par(2+il+1:2+2*il)
156 !
157 IF (present(pdx)) pdx(:) = pgrid_par(2+2*il+1:2+3*il)
158 !
159 IF (present(pdy)) pdy(:) = pgrid_par(2+3*il+1:2+4*il)
160 !
161 IF (present(kdimx)) kdimx = nint(pgrid_par(3+4*il))
162 !
163 IF (present(kdimy)) kdimy = nint(pgrid_par(4+4*il))
164 !
165 IF (present(pxall)) THEN
166  idimx= nint(pgrid_par(3+4*il))
167  pxall(:)= pgrid_par(5+4*il:4+4*il+idimx)
168 END IF
169 !
170 IF (present(pyall)) THEN
171  idimx= nint(pgrid_par(3+4*il))
172  idimy= nint(pgrid_par(4+4*il))
173  pyall(:)= pgrid_par(5+4*il+idimx:4+4*il+idimx+idimy)
174 END IF
175 !
176 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:GET_GRIDTYPE_IGN',1,zhook_handle)
177 !
178 !-------------------------------------------------------------------------------
179 END SUBROUTINE get_gridtype_ign
180 !############################################################################
181 !############################################################################
182 !############################################################################
183 ! ###################################################
184  SUBROUTINE latlon_ign(KLAMBERT,PX,PY,PLAT,PLON)
185 ! ###################################################
186 !
187 !!**** *LATLON_IGN * - Routine to compute geographical coordinates
188 !!
189 !! PURPOSE
190 !! -------
191 ! This routine computes the latitude and longitude of
192 ! an array given in LAMBERT coordinates
193 !
194 !
195 !!** METHOD
196 !! ------
197 !!
198 !! EXTERNAL
199 !! --------
200 !! None
201 !!
202 !! REFERENCE
203 !! ---------
204 !! NOTE TECHNIQUE IGN NT/G 71 :
205 !! PROJECTION CARTOGRAPHIQUE CONIQUE COMFORME DE LAMBERT
206 !! (www.ign.fr)
207 !!
208 !! AUTHOR
209 !! ------
210 !! E. Martin *Meteo-France*
211 !!
212 !! MODIFICATION
213 !! ------------
214 !! Original 10/2007
215 ! S.Lafont remplace exposant reel 2.0 par exposant entier 2
216 !-------------------------------------------------------------------------------
217 !
218 !* 0. DECLARATIONS
219 ! ------------
220 !
221 USE modd_csts,ONLY : xpi
222 USE modd_ign
223 !
224 IMPLICIT NONE
225 !
226 !* 0.1 Declarations of arguments and results
227 !
228 INTEGER, INTENT(IN) :: klambert ! Lambert type
229 REAL, DIMENSION(:), INTENT(IN) :: px,py
230  ! given conformal coordinates of the
231  ! processed points (meters);
232 REAL, DIMENSION(:), INTENT(OUT):: plat,plon
233  ! returned geographic latitudes and
234  ! longitudes of the processed points
235  ! (degrees).
236 !
237 !* 0.2 Declarations of local variables
238 REAL, DIMENSION(SIZE(PX)) :: zgamma
239 REAL, DIMENSION(SIZE(PX)) :: zr ! length of arc meridian line projection
240 REAL, DIMENSION(SIZE(PX)) :: zlatiso ! Isometric latitude
241 REAL :: zlat0 ! For iteration
242 !
243 INTEGER :: j, jj
244 REAL(KIND=JPRB) :: zhook_handle
245 !
246 !
247 !* 1. PRELIMINARY
248 ! -----------
249 !
250  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:LATLON_IGN',0,zhook_handle)
251  zr(:) =sqrt( (px(:)-xxs(klambert))**2 + (py(:)-xys(klambert))**2 )
252  zgamma(:)= atan( (px(:)-xxs(klambert)) / (xys(klambert)-py(:)) )
253 !
254 !* 2. LONGITUDE
255 ! ---------
256  plon(:)=xlonp(klambert) +zgamma(:)/xn(klambert) *180./xpi
257 !
258 !* 3. LATITUDE
259 ! --------
260  zlatiso(:)=-1./xn(klambert) * alog(abs(zr(:)/xc(klambert)))
261 !
262 !$OMP PARALLEL DO PRIVATE(JJ,J,ZLAT0)
263  DO jj=1,SIZE(plat)
264  zlat0 =2. * atan(exp(zlatiso(jj))) - xpi/2.
265  DO j=1, 1000
266  plat(jj) = 2. * atan( &
267  ( (1+xecc(klambert)*sin(zlat0))/(1-xecc(klambert)*sin(zlat0)) )**(xecc(klambert)/2.) &
268  *exp(zlatiso(jj)) ) -xpi/2.
269 !
270  IF (abs(plat(jj) - zlat0) < xcvglat ) EXIT
271  zlat0=plat(jj)
272  ENDDO
273  ENDDO
274 !$OMP END PARALLEL DO
275 !
276  plat(:)=plat(:) *180./xpi
277 !
278 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:LATLON_IGN',1,zhook_handle)
279 !---------------------------------------------------------------------------------
280 END SUBROUTINE latlon_ign
281 !---------------------------------------------------------------------------------
282 !
283 !############################################################################
284 !############################################################################
285 !############################################################################
286 !
287 ! ################################################
288  SUBROUTINE xy_ign(KLAMBERT,PX,PY,PLAT,PLON)
289 ! ################################################
290 !
291 !!**** *XY_IGN * - Routine to compute Lambert coordinates
292 !!
293 !!
294 !! PURPOSE
295 !! -------
296 ! This routine computes the Lambert coordinates
297 ! of an array given in latitude-longitude coordinates
298 !
299 !
300 !!** METHOD
301 !! ------
302 !!
303 !! WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
304 !!
305 !! EXTERNAL
306 !! --------
307 !! None
308 !!
309 !! REFERENCE
310 !! ---------
311 !! NOTE TECHNIQUE IGN NT/G 71 :
312 !! PROJECTION CARTOGRAPHIQUE CONIQUE COMFORME DE LAMBERT
313 !! (www.ign.fr)
314 !!
315 !! AUTHOR
316 !! ------
317 !! E. Martin *Meteo-France*
318 !!
319 !! MODIFICATION
320 !! ------------
321 !! Original 10/2007
322 !!
323 !-------------------------------------------------------------------------------
324 !
325 !* 0. DECLARATIONS
326 ! ------------
327 !
328 USE modd_csts, ONLY : xpi
329 USE modd_ign
330 !
331 IMPLICIT NONE
332 !
333 !* 0.1 Declarations of arguments and results
334 !
335 INTEGER, INTENT(IN) :: klambert ! Lambert type
336 REAL, DIMENSION(:), INTENT(IN) :: plat,plon
337  ! given geographic latitudes and
338  ! longitudes of the processed points
339  ! (degrees).
340 REAL, DIMENSION(:), INTENT(OUT):: px,py
341  ! returned Lambert coordinates of the
342  ! processed points (meters);
343 !
344 !* 0.2 Declarations of local variables
345 !
346 REAL :: zpi180, zpi4, zecc2
347 REAL :: zwrk ! working arrays
348 REAL :: zlatrad, zlonrad ! longitude and latitude in radians
349 REAL :: zgamma
350 REAL :: zlatfi ! Isometric latitude
351 REAL :: zr ! length of arc meridian line projection
352 INTEGER :: jj
353 REAL(KIND=JPRB) :: zhook_handle
354 !
355 !
356 !-------------------------------------------------------------------------------
357 !
358 !* 1. Latitude /Longitude in radian :
359 ! -------------------------------
360 !
361 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:XY_IGN',0,zhook_handle)
362 !
363 zpi180 = xpi / 180.
364 zpi4 = xpi / 4.
365 zecc2 = xecc(klambert) / 2.
366 !
367 !$OMP PARALLEL DO PRIVATE(JJ,ZLONRAD,ZLATRAD,ZWRK,ZLATFI,ZGAMMA,ZR)
368 DO jj=1,SIZE(plon)
369  !
370  IF (plon(jj) > 180.) THEN
371  zlonrad = (plon(jj) - 360. - xlonp(klambert)) * zpi180
372  ELSE
373  zlonrad = (plon(jj) - xlonp(klambert)) * zpi180
374  ENDIF
375  !
376  zlatrad = plat(jj) * zpi180
377  !
378 !* 2. Calcul of the isometric latitude :
379 ! ----------------------------------
380  !
381  zwrk = sin(zlatrad) * xecc(klambert)
382  !
383  zlatfi = log(tan(zpi4 + zlatrad / 2.)) + ( (log(1-zwrk)-log(1+zwrk)) * zecc2)
384  !
385 !* 3. Calcul of the lambert II coordinates X and YJJ
386 ! ---------------------------------------------
387  !
388  zr = exp(- xn(klambert) * zlatfi) * xc(klambert)
389  !
390  zgamma = xn(klambert) * zlonrad
391  !
392  px(jj) = xxs(klambert) + sin(zgamma) * zr
393  py(jj) = xys(klambert) - cos(zgamma) * zr
394  !
395 ENDDO
396 !$OMP END PARALLEL DO
397 !
398 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:XY_IGN',1,zhook_handle)
399 !-------------------------------------------------------------------------------
400 END SUBROUTINE xy_ign
401 !-------------------------------------------------------------------------------
402 !
403 !############################################################################
404 !############################################################################
405 !############################################################################
406 !
407 ! ################################################
408  SUBROUTINE map_factor_ign(KLAMBERT,PX,PY,PMAP)
409 ! ################################################
410 !
411 !!**** *MAP_FACTOR_IGN * - Routine to compute IGN map factor
412 !!
413 !!
414 !! PURPOSE
415 !! -------
416 ! Calculation of the Map factor
417 ! (ratio dist lambert / dist ellipsoide)
418 !
419 !! REFERENCE
420 !! ---------
421 !! IGN formula
422 !!
423 !! AUTHOR
424 !! ------
425 !! E. Martin *Meteo-France*
426 !!
427 !! MODIFICATION
428 !! ------------
429 !! Original 10/2007
430 ! S.Lafont remplace exposant reel 2.0 par exposant entier 2
431 !!
432 !-------------------------------------------------------------------------------
433 !
434 !* 0. DECLARATIONS
435 ! ------------
436 !
437 USE modd_csts, ONLY : xpi
438 USE modd_ign
439 !
440 IMPLICIT NONE
441 !
442 !* 0.1 Declarations of arguments and results
443 !
444 INTEGER, INTENT(IN) :: klambert ! Lambert type (1 to 6)
445 REAL, DIMENSION(:), INTENT(IN) :: px ! X lambert coordinate
446 REAL, DIMENSION(:), INTENT(IN) :: py ! Y lambert coordinate
447  ! of the processed points (m).
448 REAL, DIMENSION(:), INTENT(OUT):: pmap ! map factor
449 !
450 !* 0.2 Declarations of local variables
451 !
452 !
453 REAL, DIMENSION(SIZE(PX)) :: zlat0 ! latitude for iteration
454 REAL, DIMENSION(SIZE(PX)) :: zlat ! latitude
455 REAL, DIMENSION(SIZE(PX)) :: zlatiso ! isometric latitude
456 REAL, DIMENSION(SIZE(PX)) :: zr ! R in the IGN formula
457 REAL, DIMENSION(SIZE(PX)) :: zgrandn ! N in the IGN formula
458 !
459 INTEGER :: j
460 REAL(KIND=JPRB) :: zhook_handle
461 !
462 !-------------------------------------------------------------------------------
463 !
464 !* 1. PRELIMINARY
465 ! -----------
466 !
467  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:MAP_FACTOR_IGN',0,zhook_handle)
468  zr(:) =sqrt( (px(:)-xxs(klambert))**2 + (py(:)-xys(klambert))**2 )
469  zlatiso(:)=-1./xn(klambert) * alog(abs(zr(:)/xc(klambert)))
470  zlat0(:) =2. * atan(exp(zlatiso(:))) - xpi/2.
471 !
472  DO j=1, 100
473  zlat(:) = 2. * atan( &
474  ( (1+xecc(klambert)*sin(zlat0(:)))/(1-xecc(klambert)*sin(zlat0(:))) )**(xecc(klambert)/2.) &
475  *exp(zlatiso(:)) ) -xpi/2.
476 !
477  IF (maxval(abs(zlat(:) - zlat0(:))) < xcvglat ) EXIT
478  zlat0(:)=zlat(:)
479  ENDDO
480 !
481 !
482 !* 2. MAP FACTOR
483 ! ----------
484 !
485  zgrandn = xa(klambert) / sqrt(1-(xecc(klambert)*sin(zlat(:)))**2)
486  pmap(:)=xn(klambert)* zr(:) / ( zgrandn(:)*cos(zlat(:)) )
487 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_IGN:MAP_FACTOR_IGN',1,zhook_handle)
488 !
489 !-------------------------------------------------------------------------------
490 END SUBROUTINE map_factor_ign
491 !-------------------------------------------------------------------------------
492 !
493 
494 !############################################################################
495 !############################################################################
496 !############################################################################
497 
498 END MODULE mode_gridtype_ign
mode_gridtype_ign::xy_ign
subroutine xy_ign(KLAMBERT, PX, PY, PLAT, PLON)
Definition: mode_gridtype_ign.F90:288
modd_ign
Definition: modd_ign.F90:6
mode_gridtype_ign::get_gridtype_ign
subroutine get_gridtype_ign(PGRID_PAR, KLAMBERT, KL, PX, PY, PDX, PDY, KDIMX, KDIMY, PXALL, PYALL)
Definition: mode_gridtype_ign.F90:95
mode_gridtype_ign::latlon_ign
subroutine latlon_ign(KLAMBERT, PX, PY, PLAT, PLON)
Definition: mode_gridtype_ign.F90:184
mode_gridtype_ign
Definition: mode_gridtype_ign.F90:6
mode_gridtype_ign::map_factor_ign
subroutine map_factor_ign(KLAMBERT, PX, PY, PMAP)
Definition: mode_gridtype_ign.F90:408
mode_gridtype_ign::put_gridtype_ign
subroutine put_gridtype_ign(PGRID_PAR, KLAMBERT, PX, PY, PDX, PDY, KDIMX, KDIMY, PXALL, PYALL)
Definition: mode_gridtype_ign.F90:21
modd_csts
Definition: modd_csts.F90:6
Generated on Fri May 13 2016 09:31:04 for SURFEX by   doxygen 1.8.5