SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
mode_gridtype_conf_proj.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_conf_proj
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_conf_proj(PGRID_PAR,PLAT0,PLON0,PRPK,PBETA,&
22  plator,plonor,kimax,kjmax, &
23  px,py,pdx,pdy )
24 ! ####################################################################
25 !
26 !!**** *PUT_GRIDTYPE_CONF_PROJ* - routine to store in PGRID_PAR the horizontal grid
27 !!
28 !! AUTHOR
29 !! ------
30 !! V. Masson *Meteo France*
31 !!
32 !! MODIFICATIONS
33 !! -------------
34 !! Original 01/2004
35 !-------------------------------------------------------------------------------
36 !
37 !* 0. DECLARATIONS
38 ! ------------
39 !
40 USE modd_surf_par, ONLY : xundef
41 !
42 IMPLICIT NONE
43 !
44 !
45 !* 0.1 Declarations of arguments
46 ! -------------------------
47 !
48 REAL, INTENT(IN) :: plat0 ! reference latitude
49 REAL, INTENT(IN) :: plon0 ! reference longitude
50 REAL, INTENT(IN) :: prpk ! projection parameter
51 ! ! K=1 : stereographic north pole
52 ! ! 0<K<1 : Lambert, north hemisphere
53 ! ! K=0 : Mercator
54 ! !-1<K<0 : Lambert, south hemisphere
55 ! ! K=-1: stereographic south pole
56 REAL, INTENT(IN) :: pbeta ! angle between grid and reference longitude
57 REAL, INTENT(IN) :: plator ! latitude of point of coordinates X=0, Y=0
58 REAL, INTENT(IN) :: plonor ! longitude of point of coordinates X=0, Y=0
59 INTEGER, INTENT(IN) :: kimax ! number of points in I direction
60 INTEGER, INTENT(IN) :: kjmax ! number of points in J direction
61 REAL, DIMENSION(:), INTENT(IN) :: px ! X conformal coordinate of left boundary of grid mesh
62 REAL, DIMENSION(:), INTENT(IN) :: py ! Y conformal coordinate of bottom boundary of grid mesh
63 REAL, DIMENSION(:), INTENT(IN) :: pdx ! X grid mesh size
64 REAL, DIMENSION(:), INTENT(IN) :: pdy ! Y grid mesh size
65 REAL, DIMENSION(:), POINTER :: pgrid_par! parameters defining this grid
66 !
67 !
68 !* 0.2 Declarations of local variables
69 ! -------------------------------
70 !
71 LOGICAL :: gfull ! T : entire grid is stored
72 INTEGER :: il ! number of points
73 INTEGER :: jj ! loop counter
74 REAL(KIND=JPRB) :: zhook_handle
75 !-------------------------------------------------------------------------------
76 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:PUT_GRIDTYPE_CONF_PROJ',0,zhook_handle)
77 !
78 gfull = (kimax*kjmax==SIZE(px))
79 !
80 il = SIZE(px)
81 !
82 IF (gfull) THEN
83  !* entire grid : one can store only X and Y cooridnate arrays in each direction
84  ALLOCATE(pgrid_par(12+kimax+kjmax))
85 ELSE
86  !* only some points are present : one store the coordinates of all points
87  ALLOCATE(pgrid_par(12+2*il))
88 END IF
89 !
90 pgrid_par(1) = plat0
91 pgrid_par(2) = plon0
92 pgrid_par(3) = prpk
93 pgrid_par(4) = pbeta
94 pgrid_par(5) = plator
95 pgrid_par(6) = plonor
96 pgrid_par(7) = float(kimax)
97 pgrid_par(8) = float(kjmax)
98 IF (il>0) THEN
99  pgrid_par(9) = pdx(1)
100  pgrid_par(10)= pdy(1)
101 ELSE
102  pgrid_par(9) = xundef
103  pgrid_par(10)= xundef
104 END IF
105 pgrid_par(11) = SIZE(px)
106 IF (gfull) THEN
107  pgrid_par(12) = 1
108 ELSE
109  pgrid_par(12) = 0
110 END IF
111 !
112 IF (gfull) THEN
113  pgrid_par(12 +1:12+kimax) = px(1:kimax)
114  DO jj=1,kjmax
115  pgrid_par(12+kimax+jj) = py(1+(jj-1)*kimax)
116  END DO
117 ELSE
118  IF (il>0) THEN
119  pgrid_par(12 +1:12+ il) = px(:)
120  pgrid_par(12+ il+1:12+2*il) = py(:)
121  END IF
122 END IF
123 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:PUT_GRIDTYPE_CONF_PROJ',1,zhook_handle)
124 !-------------------------------------------------------------------------------
125 END SUBROUTINE put_gridtype_conf_proj
126 !############################################################################
127 !############################################################################
128 !############################################################################
129 ! ####################################################################
130  SUBROUTINE get_gridtype_conf_proj(PGRID_PAR,PLAT0,PLON0,PRPK,PBETA,&
131  plator,plonor,kimax,kjmax, &
132  px,py,pdx,pdy,kl )
133 ! ####################################################################
134 !
135 !!**** *GET_GRIDTYPE_CONF_PROJ* - routine to get from PGRID_PAR the horizontal grid
136 !!
137 !! AUTHOR
138 !! ------
139 !! V. Masson *Meteo France*
140 !!
141 !! MODIFICATIONS
142 !! -------------
143 !! Original 01/2004
144 !-------------------------------------------------------------------------------
145 !
146 !* 0. DECLARATIONS
147 ! ------------
148 !
149 IMPLICIT NONE
150 !
151 !
152 !* 0.1 Declarations of arguments
153 ! -------------------------
154 !
155 REAL, DIMENSION(:), INTENT(IN) :: pgrid_par! parameters defining this grid
156 REAL, INTENT(OUT), OPTIONAL :: plat0 ! reference latitude
157 REAL, INTENT(OUT), OPTIONAL :: plon0 ! reference longitude
158 REAL, INTENT(OUT), OPTIONAL :: prpk ! projection parameter
159 ! ! K=1 : stereographic north pole
160 ! ! 0<K<1 : Lambert, north hemisphere
161 ! ! K=0 : Mercator
162 ! !-1<K<0 : Lambert, south hemisphere
163 ! ! K=-1: stereographic south pole
164 REAL, INTENT(OUT), OPTIONAL :: pbeta ! angle between grid and reference longitude
165 REAL, INTENT(OUT), OPTIONAL :: plator ! latitude of point of coordinates X=0, Y=0
166 REAL, INTENT(OUT), OPTIONAL :: plonor ! longitude of point of coordinates X=0, Y=0
167 INTEGER, INTENT(OUT), OPTIONAL :: kimax ! number of points in I direction
168 INTEGER, INTENT(OUT), OPTIONAL :: kjmax ! number of points in J direction
169 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: px ! X conformal coor. of grid mesh
170 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: py ! Y conformal coor. of grid mesh
171 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: pdx ! X grid mesh size
172 REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: pdy ! Y grid mesh size
173 INTEGER, INTENT(OUT), OPTIONAL :: kl ! number of points
174 !
175 !
176 !* 0.2 Declarations of local variables
177 ! -------------------------------
178 !
179 INTEGER :: il, iimax, ijmax
180 INTEGER :: ji, jj ! loop counters
181 LOGICAL :: gfull ! T : entire grid is stored
182 !
183 REAL(KIND=JPRB) :: zhook_handle
184 !-------------------------------------------------------------------------------
185 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:GET_GRIDTYPE_CONF_PROJ',0,zhook_handle)
186 !
187 il = pgrid_par(11)
188 iimax = nint(pgrid_par(7))
189 ijmax = nint(pgrid_par(8))
190 !
191 IF (present(plat0)) plat0 = pgrid_par(1)
192 IF (present(plon0)) plon0 = pgrid_par(2)
193 IF (present(prpk )) prpk = pgrid_par(3)
194 IF (present(pbeta)) pbeta = pgrid_par(4)
195 IF (present(plator)) plator= pgrid_par(5)
196 IF (present(plonor)) plonor= pgrid_par(6)
197 IF (present(kimax)) kimax = iimax
198 IF (present(kjmax)) kjmax = ijmax
199 IF (present(pdx)) pdx(:)= pgrid_par(9)
200 IF (present(pdy)) pdy(:)= pgrid_par(10)
201 IF (present(kl)) kl = il
202 !
203 gfull = (pgrid_par(12)==1)
204 !
205 IF (present(px)) THEN
206  IF (gfull) THEN
207  DO jj=1,ijmax
208  DO ji=1,iimax
209  px(ji+(jj-1)*iimax) = pgrid_par(12+ji)
210  END DO
211  END DO
212  ELSE
213  px(:) = pgrid_par(12+1:12+il)
214  END IF
215 END IF
216 
217 IF (present(py)) THEN
218  IF (gfull) THEN
219  DO jj=1,ijmax
220  DO ji=1,iimax
221  py(ji+(jj-1)*iimax) = pgrid_par(12+iimax+jj)
222  END DO
223  END DO
224  ELSE
225  py(:) = pgrid_par(12+il+1:12+2*il)
226  END IF
227 END IF
228 !
229 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:GET_GRIDTYPE_CONF_PROJ',1,zhook_handle)
230 !
231 !-------------------------------------------------------------------------------
232 END SUBROUTINE get_gridtype_conf_proj
233 !############################################################################
234 !############################################################################
235 !############################################################################
236 ! ###################################################
237  SUBROUTINE latlon_conf_proj(PLAT0,PLON0,PRPK,PBETA,PLATOR,PLONOR, &
238  px,py,plat,plon )
239 ! ###################################################
240 !
241 !!**** *LATLON_CONF_PROJ * - Routine to compute geographical coordinates
242 !!
243 !! PURPOSE
244 !! -------
245 ! This routine computes the latitude and longitude of
246 ! an array given in cartesian conformal coordinates
247 ! Five map projections are available:
248 ! - polar-stereographic from south pole (PRPK=1),
249 ! - lambert conformal from south pole (0<PRPK<1),
250 ! - mercator (PRPK=0),
251 ! - lambert conformal from north pole (-1<PRPK<0),
252 ! - polar-stereographic from north pole (PRPK=-1).
253 !
254 !
255 !!** METHOD
256 !! ------
257 !! Spherical earth approximation is used. Longitude origin is
258 !! set in Greenwich, and is positive eastwards. An anticlockwise
259 !! rotation of PBETA degrees is applied to the conformal frame
260 !! with respect to the geographical directions.
261 !!
262 !! WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
263 !!
264 !! EXTERNAL
265 !! --------
266 !! None
267 !!
268 !! REFERENCE
269 !! ---------
270 !! Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
271 !! commun CNRM-LA, specifications techniques",
272 !! Note CNRM/GMME, 26, 139p, (Chapter 2).
273 !! Ducrocq V., 1994, "Generation de la grille dans le modele",
274 !! Note interne MNH, 5 mai, 3p.
275 !! Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
276 !! Internal note ARPEGE/ALADIN, february 27,28p.
277 !! Levallois J., 1970, "Geodesie generale", Tome 2, Collection
278 !! de l'IGN, Eyrolles, Paris, 408p.
279 !!
280 !! AUTHOR
281 !! ------
282 !! P.M. *LA*
283 !!
284 !! MODIFICATION
285 !! ------------
286 !! Original PM 24/05/94
287 !! Updated PM 27/07/94
288 !! Updated VD 23/08/94
289 !! Updated VM 24/10/95 projection from north pole (PRPK<0) and
290 !! longitudes set between PLON0-180. and PLON0+180.
291 !! Update VM 11/2004 externalized version
292 !-------------------------------------------------------------------------------
293 !
294 !* 0. DECLARATIONS
295 ! ------------
296 !
297 USE modd_csts,ONLY : xpi, xradius
298 !
299 IMPLICIT NONE
300 !
301 !* 0.1 Declarations of arguments and results
302 !
303 REAL, INTENT(IN) :: plat0 ! Reference latitude
304 REAL, INTENT(IN) :: plon0 ! Reference longitude
305 REAL, INTENT(IN) :: prpk ! projection parameter
306 ! ! K=1 : stereographic north pole
307 ! ! 0<K<1 : Lambert, north hemisphere
308 ! ! K=0 : Mercator
309 ! !-1<K<0 : Lambert, south hemisphere
310 ! ! K=-1: stereographic south pole
311 REAL, INTENT(IN) :: pbeta ! angle between grid and reference longitude
312 REAL, INTENT(IN) :: plator ! Latitude of the origine point
313 REAL, INTENT(IN) :: plonor ! Longitude of the origine point
314 REAL, DIMENSION(:), INTENT(IN) :: px,py
315  ! given conformal coordinates of the
316  ! processed points (meters);
317 REAL, DIMENSION(:), INTENT(OUT):: plat,plon
318  ! returned geographic latitudes and
319  ! longitudes of the processed points
320  ! (degrees).
321 !
322 !* 0.2 Declarations of local variables
323 !
324 REAL, DIMENSION(SIZE(PX)) :: zy
325 REAL :: zrpk,zbeta,zlat0,zlon0,zlator,zlonor
326 REAL :: zrdsdg,zclat0,zslat0,zclator,zslator
327 REAL :: zxbm0,zybm0,zro0,zga0
328 REAL :: zxp,zyp,zepsi,zt1,zcgam,zsgam,zraclat0
329 !
330 REAL, DIMENSION(SIZE(PX)) :: zata,zro2,zt2,zxmi0,zymi0
331 REAL(KIND=JPRB) :: zhook_handle
332 !
333 !--------------------------------------------------------------------------------
334 !
335 !* 1. Preliminary calculations for all projections
336 ! --------------------------------------------
337 !
338 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:LATLON_CONF_PROJ',0,zhook_handle)
339 zrdsdg = xpi/180. ! Degree to radian conversion factor
340 zepsi = 10.*epsilon(1.) ! A small number
341 !
342 ! By definition, (PLONOR,PLATOR) are the geographical
343 ! coordinates, and (ZXBM0,ZYBM0) the conformal cartesian
344 ! point of coordinates (x=0,y=0).
345 !
346 zxbm0 = 0.
347 zybm0 = 0.
348 !
349 !-------------------------------------------------------------------------------
350 !
351 !* 2. POLAR STEREOGRAPHIC AND LAMBERT CONFORMAL CASES
352 ! -----------------------------------------------
353 ! (PRPK=1 P-stereo, 0<PRPK<1 Lambert)
354 !
355 IF(prpk /= 0.) THEN
356 !
357  IF (prpk<0.) THEN ! projection from north pole
358  zrpk=-prpk
359  zbeta=-pbeta
360  zlat0=-plat0
361  zlon0=plon0+180.
362  zlator=-plator
363  zlonor=plonor+180.
364  zy(:)=-py(:)
365  zybm0=-zybm0
366  ELSE ! projection from south pole
367  zrpk=prpk
368  zbeta=pbeta
369  zlat0=plat0
370  zlon0=plon0
371  zlator=plator
372  zlonor=plonor
373  zy(:)=py(:)
374  ENDIF
375 !
376 !* 2.1 Preliminary calculations
377 !
378  zclat0 = cos(zrdsdg*zlat0)
379  zslat0 = sin(zrdsdg*zlat0)
380  zclator = cos(zrdsdg*zlator)
381  zslator = sin(zrdsdg*zlator)
382  zro0 = (xradius/zrpk)*(abs(zclat0))**(1.-zrpk) &
383  * ((1.+zslat0)*abs(zclator)/(1.+zslator))**zrpk
384  zga0 = (zrpk*(zlonor-zlon0)-zbeta)*zrdsdg
385  zxp = zxbm0-zro0*sin(zga0)
386  zyp = zybm0+zro0*cos(zga0)
387 !
388 !* 2.2 Longitude
389 !
390  WHERE (abs(zy(:)-zyp) < zepsi &
391  .AND.abs(px(:)-zxp) < zepsi)
392  zata(:) = 0.
393  ELSEWHERE
394  zata(:) = atan2(-(zxp-px(:)),(zyp-zy(:)))/zrdsdg
395  END WHERE
396  !
397  plon(:) = (zbeta+zata(:))/zrpk+zlon0
398 !
399 !* 2.3 Latitude
400 !
401  zro2(:) = (px(:)-zxp)**2+(zy(:)-zyp)**2
402  zt1 = (xradius*(abs(zclat0))**(1.-zrpk))**(2./zrpk) &
403  * (1+zslat0)**2
404  zt2(:) = (zrpk**2*zro2(:))**(1./zrpk)
405  !
406  plat(:) = (xpi/2.-acos((zt1-zt2(:))/(zt1+zt2(:))))/zrdsdg
407 !
408  IF (prpk<0.) THEN ! projection from north pole
409  plat(:)=-plat(:)
410  plon(:)=plon(:)+180.
411  ENDIF
412 !
413 !-------------------------------------------------------------------------------
414 !
415 !* 3. MERCATOR PROJECTION WITH ROTATION
416 ! ---------------------------------
417 ! (PRPK=0)
418 !
419 ELSE
420 !
421 !* 3.1 Preliminary calculations
422 !
423  zcgam = cos(-zrdsdg*pbeta)
424  zsgam = sin(-zrdsdg*pbeta)
425  zraclat0 = xradius*cos(zrdsdg*plat0)
426 !
427 !* 3.2 Longitude
428 !
429  zxmi0(:) = px(:)-zxbm0
430  zymi0(:) = py(:)-zybm0
431  !
432  plon(:) = (zxmi0(:)*zcgam+zymi0(:)*zsgam) &
433  / (zraclat0*zrdsdg)+plonor
434 !
435 !* 3.3 Latitude
436 !
437  zt1 = alog(tan(xpi/4.+plator*zrdsdg/2.))
438  zt2(:) = (-zxmi0(:)*zsgam+zymi0(:)*zcgam)/zraclat0
439  !
440  plat(:) = (-xpi/2.+2.*atan(exp(zt1+zt2(:))))/zrdsdg
441 !
442 !---------------------------------------------------------------------------------
443 !
444 !* 4. EXIT
445 ! ----
446 !
447 END IF
448 plon(:)=plon(:)+nint((plon0-plon(:))/360.)*360.
449 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:LATLON_CONF_PROJ',1,zhook_handle)
450 !---------------------------------------------------------------------------------
451 END SUBROUTINE latlon_conf_proj
452 !---------------------------------------------------------------------------------
453 !
454 !############################################################################
455 !############################################################################
456 !############################################################################
457 !
458 ! ################################################
459  SUBROUTINE xy_conf_proj(PLAT0,PLON0,PRPK,PBETA,PLATOR,PLONOR, &
460  px,py,plat,plon )
461 ! ################################################
462 !
463 !!**** *XY_CONF_PROJ * - Routine to compute conformal coordinates
464 !!
465 !!
466 !! PURPOSE
467 !! -------
468 ! This routine computes the cartesian conformal coordinates
469 ! of an array given in latitude-longitude coordinates
470 ! Three map projections are available:
471 ! - polar-stereographic (XRPK=1),
472 ! - lambert conformal (0<XRPK<1),
473 ! - mercator (XRPK=0).
474 !
475 !
476 !!** METHOD
477 !! ------
478 !! Spherical earth approximation is used. Longitude origin is
479 !! set in Greenwich, and is positive eastwards. An anticlockwise
480 !! rotation of XBETA degrees is applied to the conformal frame
481 !! with respect to the geographical directions.
482 !!
483 !! WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
484 !!
485 !! EXTERNAL
486 !! --------
487 !! None
488 !!
489 !! REFERENCE
490 !! ---------
491 !! Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
492 !! commun CNRM-LA, specifications techniques",
493 !! Note CNRM/GMME, 26, 139p, (Chapter 2).
494 !! Ducrocq V., 1994, "Generation de la grille dans le modele",
495 !! Note interne MNH, 5 mai, 3p.
496 !! Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
497 !! Internal note ARPEGE/ALADIN, february 27,28p.
498 !! Levallois J., 1970, "Geodesie generale", Tome 2, Collection
499 !! de l'IGN, Eyrolles, Paris, 408p.
500 !!
501 !! AUTHOR
502 !! ------
503 !! P.M. *LA*
504 !!
505 !! MODIFICATION
506 !! ------------
507 !! Original PM 24/05/94
508 !! Updated PM 27/07/94
509 !! Updated VD 23/08/94
510 !! Updated VM 24/10/95 projection from north pole (XRPK<0) and
511 !! longitudes set between XLON0-180. and XLON0+180.
512 !! Updated VM 03/2004 externalized version
513 !!
514 !-------------------------------------------------------------------------------
515 !
516 !* 0. DECLARATIONS
517 ! ------------
518 !
519 USE modd_csts, ONLY : xpi, xradius
520 !
521 IMPLICIT NONE
522 !
523 !* 0.1 Declarations of arguments and results
524 !
525 REAL, INTENT(IN) :: plat0 ! Reference latitude
526 REAL, INTENT(IN) :: plon0 ! Reference longitude
527 REAL, INTENT(IN) :: prpk ! projection parameter
528 ! ! K=1 : stereographic north pole
529 ! ! 0<K<1 : Lambert, north hemisphere
530 ! ! K=0 : Mercator
531 ! !-1<K<0 : Lambert, south hemisphere
532 ! ! K=-1: stereographic south pole
533 REAL, INTENT(IN) :: pbeta ! angle between grid and reference longitude
534 REAL, INTENT(IN) :: plator ! Latitude of the origine point
535 REAL, INTENT(IN) :: plonor ! Longitude of the origine point
536 REAL, DIMENSION(:), INTENT(OUT):: px,py
537  ! returned conformal coordinates of the
538  ! processed points (meters);
539 REAL, DIMENSION(:), INTENT(IN) :: plat,plon
540  ! given geographic latitudes and
541  ! longitudes of the processed points
542  ! (degrees).
543 !
544 !* 0.2 Declarations of local variables
545 !
546 REAL,DIMENSION(SIZE(PLAT)) :: zlat,zlon
547 REAL :: zrpk,zbeta,zlat0,zlon0,zlator,zlonor
548 REAL :: zrdsdg,zclat0,zslat0,zclator,zslator
549 REAL :: zxbm0,zybm0,zro0,zga0
550 REAL :: zxp,zyp,zcgam,zsgam,zraclat0,zxe,zye
551 !
552 REAL,DIMENSION(SIZE(PLAT)) :: zclat,zslat,zro,zga,zxpr,zypr
553 REAL(KIND=JPRB) :: zhook_handle
554 !
555 !
556 !-------------------------------------------------------------------------------
557 !
558 !* 1. PRELIMINARY CALCULATION FOR ALL PROJECTIONS
559 ! -------------------------------------------
560 !
561 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:XY_CONF_PROJ',0,zhook_handle)
562 zrdsdg = xpi/180. ! Degree to radian conversion factor
563 !
564 ! By definition, (PLONOR,PLATOR) are the geographical
565 ! coordinates, and (ZXBM0,ZYBM0) the conformal cartesian
566 ! coordinates of the x=0, y=0 point.
567 !
568 zxbm0 = 0.
569 zybm0 = 0.
570 !
571 zlon(:)=plon(:)
572 zlon(:)=zlon(:)+nint((plon0-zlon(:))/360.)*360.
573 !
574 zlonor=plonor
575 zlonor=zlonor+nint((plon0-zlonor)/360.)*360.
576 !------------------------------------------------------------------------------
577 !
578 !* 2. POLAR SEREOGRAPHIC AND LAMBERT CONFORMAL CASES
579 ! ----------------------------------------------
580 ! (PRPK=1 P-stereo, 0<PRPK<1 Lambert)
581 !
582 IF(prpk /= 0.) THEN
583 !
584  IF (prpk<0.) THEN ! projection from north pole
585  zrpk=-prpk
586  zbeta=-pbeta
587  zlat0=-plat0
588  zlon0=plon0+180.
589  zlator=-plator
590  zlonor=zlonor+180.
591  zlat(:)=-plat(:)
592  zlon(:)=zlon(:)+180.
593  zybm0=-zybm0
594  ELSE ! projection from south pole
595  zrpk=prpk
596  zbeta=pbeta
597  zlat0=plat0
598  zlon0=plon0
599  zlator=plator
600  zlonor=zlonor
601  zlat(:)=plat(:)
602  zlon(:)=zlon(:)
603  ENDIF
604 !
605 !* 2.1 Preliminary calculations
606 !
607  zclat0 = cos(zrdsdg*zlat0)
608  zslat0 = sin(zrdsdg*zlat0)
609  zclator = cos(zrdsdg*zlator)
610  zslator = sin(zrdsdg*zlator)
611  zro0 = (xradius/zrpk)*(abs(zclat0))**(1.-zrpk) &
612  * ((1.+zslat0)*abs(zclator)/(1.+zslator))**zrpk
613  zga0 = (zrpk*(zlonor-zlon0)-zbeta)*zrdsdg
614  zxp = zxbm0-zro0*sin(zga0)
615  zyp = zybm0+zro0*cos(zga0)
616 !
617 !* 2.2 Conformal coordinates in meters
618 !
619  zclat(:) = cos(zrdsdg*zlat(:))
620  zslat(:) = sin(zrdsdg*zlat(:))
621  zro(:) = (xradius/zrpk)*(abs(zclat0))**(1.-zrpk) &
622  * ((1.+zslat0)*abs(zclat(:))/(1.+zslat(:)))**zrpk
623  zga(:) = (zrpk*(zlon(:)-zlon0)-zbeta)*zrdsdg
624 !
625  px(:) = zxp+zro(:)*sin(zga(:))
626  py(:) = zyp-zro(:)*cos(zga(:))
627 !
628  IF (prpk<0.) THEN ! projection from north pole
629  py(:)=-py(:)
630  ENDIF
631 !
632 !-------------------------------------------------------------------------------
633 !
634 !* 3. MERCATOR PROJECTION WITH ROTATION
635 ! ---------------------------------
636 ! (PRPK=0)
637 !
638 ELSE
639 !
640 !* 3.1 Preliminary calculations
641 !
642  zcgam = cos(-zrdsdg*pbeta)
643  zsgam = sin(-zrdsdg*pbeta)
644  zraclat0 = xradius*cos(zrdsdg*plat0)
645  zxe = zxbm0*zcgam+zybm0*zsgam &
646  - zraclat0*(plonor-plon0)*zrdsdg
647  zye =-zxbm0*zsgam+zybm0*zcgam &
648  - zraclat0*log(tan(xpi/4.+plator*zrdsdg/2.))
649 !
650 !* 3.2 Conformal coordinates
651 !
652  zxpr(:) = zraclat0*(zlon(:)-plon0)*zrdsdg+zxe
653  zypr(:) = zraclat0*log(tan(xpi/4.+plat(:)*zrdsdg/2.))+zye
654  !
655  px(:) = zxpr(:)*zcgam-zypr(:)*zsgam
656  py(:) = zxpr(:)*zsgam+zypr(:)*zcgam
657 !
658 !-------------------------------------------------------------------------------
659 !
660 !* 4. EXIT
661 ! ----
662 !
663 END IF
664 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:XY_CONF_PROJ',1,zhook_handle)
665 !-------------------------------------------------------------------------------
666 END SUBROUTINE xy_conf_proj
667 !-------------------------------------------------------------------------------
668 !
669 !############################################################################
670 !############################################################################
671 !############################################################################
672 !
673 ! ################################################
674  SUBROUTINE map_factor_conf_proj(PLAT0,PRPK,PLAT,PMAP)
675 ! ################################################
676 !
677 !!**** *MAP_FACTOR_CONF_PROJ * - Routine to compute conformal coordinates
678 !!
679 !!
680 !! PURPOSE
681 !! -------
682 !
683 !! REFERENCE
684 !! ---------
685 !! Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
686 !! commun CNRM-LA, specifications techniques",
687 !! Note CNRM/GMME, 26, 139p, (Chapter 2).
688 !! Ducrocq V., 1994, "Generation de la grille dans le modele",
689 !! Note interne MNH, 5 mai, 3p.
690 !! Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
691 !! Internal note ARPEGE/ALADIN, february 27,28p.
692 !! Levallois J., 1970, "Geodesie generale", Tome 2, Collection
693 !! de l'IGN, Eyrolles, Paris, 408p.
694 !!
695 !! AUTHOR
696 !! ------
697 !! P.M. *LA*
698 !!
699 !! MODIFICATION
700 !! ------------
701 !! Original PM 24/05/94
702 !! Updated PM 27/07/94
703 !! Updated VD 23/08/94
704 !! Updated VM 24/10/95 projection from north pole (XRPK<0) and
705 !! longitudes set between XLON0-180. and XLON0+180.
706 !! Updated VM 03/2004 externalized version
707 !!
708 !-------------------------------------------------------------------------------
709 !
710 !* 0. DECLARATIONS
711 ! ------------
712 !
713 USE modd_csts, ONLY : xpi, xradius
714 !
715 IMPLICIT NONE
716 !
717 !* 0.1 Declarations of arguments and results
718 !
719 REAL, INTENT(IN) :: plat0 ! Reference latitude
720 REAL, INTENT(IN) :: prpk ! projection parameter
721 ! ! K=1 : stereographic north pole
722 ! ! 0<K<1 : Lambert, north hemisphere
723 ! ! K=0 : Mercator
724 ! !-1<K<0 : Lambert, south hemisphere
725 ! ! K=-1: stereographic south pole
726 REAL, DIMENSION(:), INTENT(IN) :: plat ! given geographic latitudes
727  ! of the processed points (degrees).
728 REAL, DIMENSION(:), INTENT(OUT):: pmap ! map factor
729 !
730 !* 0.2 Declarations of local variables
731 !
732 !
733 REAL :: zlat0 ! reference latitude
734 REAL :: zrpk ! projection parameter
735 ! ! K=1 : stereographic north pole
736 ! ! 0<K<1 : Lambert, north hemisphere
737 ! ! K=0 : Mercator
738 ! !-1<K<0 : Lambert, south hemisphere
739 ! ! K=-1: stereographic south pole
740 REAL, DIMENSION(SIZE(PLAT)) :: zlat ! latitude
741 !
742 REAL :: zclat0 ! cos(lat0)
743 REAL :: zslat0 ! sin(lat0)
744 REAL :: zrdsdg ! pi/180
745 LOGICAL :: gnorthproj! T: projection from north pole
746 REAL(KIND=JPRB) :: zhook_handle
747 !
748 !-------------------------------------------------------------------------------
749 !
750 !* 1. PRELIMINARY CALCULATION FOR ALL PROJECTIONS
751 ! -------------------------------------------
752 !
753 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:MAP_FACTOR_CONF_PROJ',0,zhook_handle)
754 zrdsdg = xpi/180. ! Degree to radian conversion factor
755 !
756 gnorthproj = prpk < 0.
757 IF (gnorthproj) THEN ! projection from north pole
758  zrpk=-prpk
759  zlat0=-plat0
760  zlat(:)=-plat(:)
761 ELSE
762  zrpk=prpk
763  zlat0=plat0
764  zlat(:)=plat(:)
765 ENDIF
766 !
767 zrdsdg = xpi/180.
768 zclat0 = cos(zrdsdg*zlat0)
769 zslat0 = sin(zrdsdg*zlat0)
770 !
771 IF (abs(zclat0)<1.e-10 .AND. (abs(zrpk-1.)<1.e-10)) THEN
772  pmap(:) = (1.+zslat0)/(1.+sin(zrdsdg*zlat(:)))
773 ELSE
774  WHERE (abs(cos(zrdsdg*zlat(:)))>1.e-10)
775  pmap(:) = ((zclat0/cos(zrdsdg*zlat(:)))**(1.-zrpk)) &
776  * ((1.+zslat0)/(1.+sin(zrdsdg*zlat(:))))**zrpk
777  ELSEWHERE
778  pmap(:) = (1.+zslat0)/(1.+sin(zrdsdg*zlat(:)))
779  ENDWHERE
780 END IF
781 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_CONF_PROJ:MAP_FACTOR_CONF_PROJ',1,zhook_handle)
782 !
783 !-------------------------------------------------------------------------------
784 END SUBROUTINE map_factor_conf_proj
785 !-------------------------------------------------------------------------------
786 !
787 !############################################################################
788 !############################################################################
789 !############################################################################
790 
791 END MODULE mode_gridtype_conf_proj
mode_gridtype_conf_proj::map_factor_conf_proj
subroutine map_factor_conf_proj(PLAT0, PRPK, PLAT, PMAP)
Definition: mode_gridtype_conf_proj.F90:674
mode_gridtype_conf_proj::put_gridtype_conf_proj
subroutine put_gridtype_conf_proj(PGRID_PAR, PLAT0, PLON0, PRPK, PBETA, PLATOR, PLONOR, KIMAX, KJMAX, PX, PY, PDX, PDY)
Definition: mode_gridtype_conf_proj.F90:21
mode_gridtype_conf_proj::latlon_conf_proj
subroutine latlon_conf_proj(PLAT0, PLON0, PRPK, PBETA, PLATOR, PLONOR, PX, PY, PLAT, PLON)
Definition: mode_gridtype_conf_proj.F90:237
mode_gridtype_conf_proj::xy_conf_proj
subroutine xy_conf_proj(PLAT0, PLON0, PRPK, PBETA, PLATOR, PLONOR, PX, PY, PLAT, PLON)
Definition: mode_gridtype_conf_proj.F90:459
modd_surf_par
Definition: modd_surf_par.F90:6
mode_gridtype_conf_proj::get_gridtype_conf_proj
subroutine get_gridtype_conf_proj(PGRID_PAR, PLAT0, PLON0, PRPK, PBETA, PLATOR, PLONOR, KIMAX, KJMAX, PX, PY, PDX, PDY, KL)
Definition: mode_gridtype_conf_proj.F90:130
mode_gridtype_conf_proj
Definition: mode_gridtype_conf_proj.F90:6
modd_csts
Definition: modd_csts.F90:6
Generated on Fri May 13 2016 09:31:04 for SURFEX by   doxygen 1.8.5