Philippe 02/05/2016: moved all LIBTOOLS files in LIBTOOLS directory
[MNH-git_open_source-lfs.git] / LIBTOOLS / tools / diachro / src / mesonh_MOD / mode_gridproj.f90
1 !     ####################
2       MODULE MODE_GRIDPROJ
3 !     ####################
4 !
5 !!****  *MODE_GRIDPROJ*  -   module routine  SM_GRIDPROJ
6 !!
7 !!      PURPOSE
8 !!      -------
9 !         This executable module packages a set of cartographic
10 !       module-procedures: 
11 !
12 !       SM_GRIDPROJ : to  compute the Jacobian in the case of 
13 !                    conformal projection;
14 !       SM_LATLON   : to compute geographic  from conformal
15 !                    cartesian coordinates;
16 !       SM_XYHAT    : to compute conformal cartesian from
17 !                     geographic coordinates;
18 !       LATREF2     : to compute the second reference latitude
19 !                    in the case of Lambert conformal projection
20 !
21 !!
22 !!**    IMPLICIT ARGUMENTS
23 !!      ------------------
24 !!           NONE
25 !!
26 !!      AUTHOR
27 !!      ------
28 !!          P.M.         *LA*
29 !!
30 !!      MODIFICATION
31 !!      ------------
32 !!          Original  24/05/94
33 !!
34 !!    
35 !------------------------------------------------------------------------------
36 !
37 !*                0.  DECLARATIONS
38 !                     ------------
39 !------------------------------------------------------------------------------
40 !
41 INTERFACE SM_LATLON
42    MODULE PROCEDURE SM_LATLON_A,SM_LATLON_S
43 END INTERFACE
44 INTERFACE SM_XYHAT
45    MODULE PROCEDURE SM_XYHAT_A,SM_XYHAT_S
46 END INTERFACE
47 !
48 CONTAINS
49 !-------------------------------------------------------------------------------
50 !-------------------------------------------------------------------------------
51 !
52 !*                1.  ROUTINE  SM_GRIDPROJ   
53 !                     --------------------
54 !-------------------------------------------------------------------------------
55 !      ####################################################################
56        SUBROUTINE SM_GRIDPROJ(HLUOUT,PXHAT,PYHAT,PZHAT,PZS,           &
57                               OSLEVE,PLEN1,PLEN2,PZSMT,PLATOR,PLONOR, &
58                               PMAP,PLAT,PLON,PDXHAT,PDYHAT,PZZ,PJ)
59 !      ####################################################################
60 !
61 !!*****  *SM_GRIDPROJ * - Computes Jacobian J, map factor M,
62 !!    horizontal grid-meshes, latitude and longitude  at the 
63 !!    "mass" point locations.
64 !!
65 !!
66 !!    PURPOSE
67 !!    -------
68 !       The purpose of this routine is to compute the Jacobian (J) at the
69 !     "mass" point location in the case of a conformal projection.
70 !     The map factor of the projection, the horizontal mesh-sizezs, and the 
71 !     the geograpical locations are also computed in the course of 
72 !     this calculation.
73 !        Five map projections are available: 
74 !      - polar-stereographic from south pole  (XRPK=1),
75 !      - lambert conformal from south pole  (0<XRPK<1),
76 !      - mercator                             (XRPK=0),
77 !      - lambert conformal from north pole (-1<XRPK<0),
78 !      - polar-stereographic from north pole  (XRPK=-1).
79 !
80 !
81 !!**   METHOD
82 !!     ------
83 !!       The height, and the correction for spherical earth are first computed.
84 !!     Next, the conformal horizontal locations, the geographical coordinates 
85 !!     and the map factor  are derived at the "mass" grid-points.  
86 !!     The same formula can (hopefully) be used for all the projections cases
87 !!     (see Joly, 1992).
88 !!
89 !!       WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
90 !!
91 !!
92 !!    EXTERNAL
93 !!    --------
94 !!      NONE
95 !!
96 !!    EXPLICIT ARGUMENTS (not required, but given for convenience)
97 !!    ------------------
98 !!       PXHAT   : conformal coordinate x  (meters, u-grid, input)
99 !!       PYHAT   : conformal coordinate y  (meters, v-grid, input)
100 !!       PZHAT   : Gal-chen altitude  zhat (meters, w-grid, input)
101 !!       PZS     : topography              (meters, masss-grid, input)
102 !!       PLATOR  : Latitude of the origine point (degrees, mass grid, input)
103 !!       PLONOR  : Longitude of the origine point (degrees, mass grid, input)
104 !!       PMAP    : map scale               (no-unit, mass-grid, output)
105 !!       PLAT    : latitude                (degrees, mass-grid, output)
106 !!       PLON    : longitude               (degrees, mass-grid, output)
107 !!       PDXHAT  : local x mesh size       (meters, u-grid, output)
108 !!       PDYHAT  : local y mesh size       (meters, v-grid, output)
109 !!       PZZ     : true altitude z         (meters, w-grid, output)
110 !!       PJ      : jacobian                (no-unit, mass-grid, output)
111 !!
112 !!    IMPLICIT ARGUMENTS
113 !!    ------------------
114 !!       Module MODD_CONF      : contains declaration of configuration variables 
115 !!          LTHINSHELL   : Logical for thinshell approximation
116 !!          NVERB        : Listing verbosity control
117 !!
118 !!       Module MODD_CST       : contains Physical constants
119 !!          XPI          : Pi;    
120 !!          XRADIUS      : Earth radius (meters);
121 !!
122 !!       Module MODD_PARAMETERS: contains declaration of parameter variables 
123 !!          JPHEXT       : horizontal depth of arrays border 
124 !!          JPVEXT       : vertical   depth of arrays border
125 !!
126 !!       Module MODD_GRID      : contains spatial grid variables
127 !!          XLAT0   : map reference latitude  (degrees)
128 !!          XRPK    : projection parameter    (no-unit)
129 !!                            
130 !!
131 !!    REFERENCE
132 !!    ---------
133 !!      Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
134 !!            commun CNRM-LA, specifications techniques", 
135 !!            Note CNRM/GMME, 26, 139p, (Chapter 2).
136 !!      Ducrocq V., 1994, "Generation de la grille dans le modele",
137 !!            Note interne MNH, 5 mai, 3p.
138 !!      Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
139 !!            Internal note ARPEGE/ALADIN, february 27,28p.
140 !!      Levallois J., 1970, "Geodesie generale", Tome 2, Collection
141 !!             de l'IGN, Eyrolles, Paris, 408p.
142 !!      (chapters 2 and 3)
143 !!
144 !!
145 !!    AUTHOR
146 !!    ------
147 !!      P. Mascart        * LA *
148 !!
149 !!    MODIFICATIONS
150 !!    -------------
151 !!      Original    PM  20/06/94  (from SM_GRIDCART by V. Ducrocq)
152 !!      Updated     PM  26/07/94
153 !!      Updated     VD  23/08/94
154 !!                      14/04/95  (Masson) bug in the ZYHTAM computation 
155 !!                      24/10/95  (Masson) controls during PMAP computation and
156 !!                                         projection from north pole (XPRK<0)
157 !!                      14/03/96  (Masson) enforce  -180<XLONOR<+180     
158 !!                      01/11/96  (Mallet) bug for the MAP FACTOR computation
159 !!      Sleve coordinate        G. Zangler  *LA*             nov 2005
160 !!!
161 !-------------------------------------------------------------------------------
162 !
163 !*       0.    DECLARATIONS
164 !              ------------
165 !
166 USE MODD_CONF          
167 USE MODD_CST          
168 USE MODD_PARAMETERS 
169 USE MODD_GRID      
170 !
171 USE MODI_VERT_COORD
172 !
173 IMPLICIT NONE
174 !
175 !*       0.1   Declarations of arguments
176 !
177 CHARACTER(LEN=*),        INTENT(IN) :: HLUOUT            ! Name of output-listing
178 REAL, DIMENSION(:),      INTENT(IN) :: PXHAT,PYHAT,PZHAT ! Positions x,y,z in 
179                                                          ! the cartesian plane
180 REAL, DIMENSION(:,:),    INTENT(IN) :: PZS               ! Orography
181 LOGICAL,                INTENT(IN)  :: OSLEVE          ! flag for SLEVE coordinate
182 REAL,                   INTENT(IN)  :: PLEN1           ! Decay scale for smooth topography
183 REAL,                   INTENT(IN)  :: PLEN2           ! Decay scale for small-scale topography deviation
184 REAL, DIMENSION(:,:),   INTENT(IN)  :: PZSMT           ! smooth orography
185 REAL,                    INTENT(IN) :: PLATOR            ! Latitude of the 
186                                                          ! origine point 
187 REAL,                    INTENT(IN) :: PLONOR            ! Longitude of the 
188                                                          ! origine point 
189 REAL, DIMENSION(:),     INTENT(OUT) :: PDXHAT          ! Local meshlength in 
190                                                        ! x direction
191 REAL, DIMENSION(:),     INTENT(OUT) :: PDYHAT          ! Local meshlength in 
192                                                        ! y direction 
193 REAL, DIMENSION(:,:),   INTENT(OUT) :: PMAP            ! Local map scale
194                                                        ! of mass gridpoints
195 REAL, DIMENSION(:,:),   INTENT(OUT) :: PLAT,PLON       ! Latitude and longitude
196                                                        ! of mass gridpoints
197 REAL, DIMENSION(:,:,:), INTENT(OUT):: PZZ              ! True altitude of the
198                                                        ! w grid-point
199 REAL, DIMENSION(:,:,:), INTENT(OUT):: PJ               ! Jacobian of the
200                                                        ! GCS transformation
201                                                        ! of mass gridpoints
202 !
203 !*       0.2   Declarations of local variables
204 !
205 REAL, DIMENSION(SIZE(PXHAT,1),SIZE(PYHAT,1),SIZE(PZHAT,1)):: ZDZ ! Local z
206                                                                  ! meshsize
207 REAL                                        :: ZH       ! H 
208 REAL, DIMENSION(SIZE(PXHAT,1),SIZE(PYHAT,1)):: ZCOEF    ! 1-zs/H 
209                                                         ! upper bounds in 
210 REAL, DIMENSION(SIZE(PXHAT,1),SIZE(PYHAT,1),SIZE(PZHAT,1)):: ZAPZOA2 ! Spherical
211                                                                      ! earth factor
212                                                                      ! for J  
213 REAL, DIMENSION(SIZE(PXHAT,1),SIZE(PYHAT,1)):: ZXHATM   ! X and Y mass point
214 REAL, DIMENSION(SIZE(PXHAT,1),SIZE(PYHAT,1)):: ZYHATM   ! conformal coordinates
215
216 REAL ZRDSDG                              ! Radian to Degree conversion factor
217 REAL ZCLAT0,ZSLAT0                       ! Cos and Sin of XLAT0
218 REAL,DIMENSION(SIZE(PLAT,1),SIZE(PLAT,2)) :: ZLAT
219 REAL                                      :: ZRPK,ZLAT0
220 !
221 INTEGER      :: IIU,IJU,IKU      ! Uupper bounds of PXHAT,PYHAT,PZHAT  
222 INTEGER      :: IIE,IJE,IKE      ! End of usefull area of PXHAT,PYHAT,PZHAT  
223 INTEGER      :: IIB,IJB,IKB      ! Begining of usefull area of PXHAT,PYHAT,PZHAT  
224 INTEGER      :: IDELTA1          ! Switch=0 if thin shell approximation
225 INTEGER      :: ILUOUT,IRESP     ! Unit number for prints, FM error code 
226 INTEGER      :: JKLOOP           ! Index for control prints
227 !
228 !-------------------------------------------------------------------------------
229 !
230 !*       1.    RETRIEVE LOGICAL UNIT NUMBER FOR OUTPUT-LISTING AND DIMENSIONS 
231 !              --------------------------------------------------------------
232 !
233 CALL FMLOOK(HLUOUT,HLUOUT,ILUOUT,IRESP)
234 !
235 IIU = UBOUND(PXHAT,1)     
236 IJU = UBOUND(PYHAT,1)    
237 IKU = UBOUND(PZHAT,1)       
238 IIE = IIU-JPHEXT
239 IJE = IJU-JPHEXT
240 IKE = IKU-JPVEXT
241 IIB = 1+JPHEXT
242 IJB = 1+JPHEXT
243 IKB = 1+JPVEXT
244 !
245 IF(NVERB >= 10) THEN                               !Value control
246   WRITE(ILUOUT,*) 'SM_GRIDPROJ: IIU,IJU,IKU=',IIU,IJU,IKU
247   WRITE(ILUOUT,*) 'SM_GRIDPROJ: IIE,IJE,IKE=',IIE,IJE,IKE
248   WRITE(ILUOUT,*) 'SM_GRIDPROJ: IIB,IJB,IKB=',IIB,IJB,IKB
249 END IF
250 !
251 !-------------------------------------------------------------------------------
252 !
253 !*       2.    COMPUTES Z   (W LEVEL)
254 !              ----------------------
255 !
256 !JDJDJDJD 291196
257 ! Ai enleve le forcage ci-apres --> non compatibilite avec la partie CONVERSION
258 ! actuelle
259 !CSTORAGE_TYPE='PG'
260 !print *,' MODE_GRIDPROJ CSTORAGE_TYPE AP FORCAGE TEMPORAIRE ',CSTORAGE_TYPE
261 IF((CCONF /= 'POSTP') .OR. (CCONF =='POSTP' .AND. CSTORAGE_TYPE /= 'PG' &
262                                             .AND. CSTORAGE_TYPE /= 'SU' ))THEN
263 !JDJDJDJD 291196
264 !
265 CALL VERT_COORD(OSLEVE,PZS,PZSMT,PLEN1,PLEN2,PZHAT,PZZ)
266 !
267 IF(NVERB >= 10) THEN                               !Value control
268   WRITE(ILUOUT,*) 'SM_GRIDPROJ: Some PZS values:'
269   WRITE(ILUOUT,*) PZS(1,1),PZS(IIU/2,IJU/2),PZS(IIU,IJU)  
270   WRITE(ILUOUT,*) 'SM_GRIDPROJ: Some PZZ values:'
271   DO JKLOOP=1,IKU
272     WRITE(ILUOUT,*) PZZ(1,1,JKLOOP),PZZ(IIU/2,IJU/2,JKLOOP), &
273                     PZZ(IIU,IJU,JKLOOP)  
274   END DO
275 END IF
276 !
277 !-------------------------------------------------------------------------------
278 !
279 !*       3.   COMPUTE SPHERICAL EARTH FACTOR (MASS LEVEL)
280 !             --------------------------------------------
281 !
282 !     NOTE: In this routine LCARTESIAN is ALWAYS .F.
283 !           Hence, IDELTA2 is always set to 1
284 !
285 IF     (LTHINSHELL) IDELTA1=0                ! THIN SHELL APPROX.
286 IF(.NOT.LTHINSHELL) IDELTA1=1                ! NO THIN SHELL APPROX.
287 !
288 IF(NVERB >= 10) THEN                         !Value control
289   WRITE(ILUOUT,*) 'SM_GRIDPROJ: LTHINSHELL, IDELTA1=',LTHINSHELL,IDELTA1
290   WRITE(ILUOUT,*) 'SM_GRIDPROJ: XRADIUS=',XRADIUS
291 ENDIF
292 !
293 ! For the time being, an inline implementation of  MZF
294 ! is provided here.
295 !
296 ZAPZOA2(:,:,1:IKU-1) = (.5*((XRADIUS+IDELTA1*PZZ(:,:,1:IKU-1))    &
297                      + (XRADIUS+IDELTA1*PZZ(:,:,2:IKU)))          &
298                        /XRADIUS)**2
299 ZAPZOA2(:,:,IKU)     = 2.*ZAPZOA2(:,:,IKU-1)-ZAPZOA2(:,:,IKU-2)
300 !
301 IF(NVERB >= 10) THEN                               !Value control
302   WRITE(ILUOUT,*) 'SM_GRIDPROJ: Some ZAPZOA2 values:'
303   DO JKLOOP=1,IKU
304     WRITE(ILUOUT,*) ZAPZOA2(1,1,JKLOOP),ZAPZOA2(IIU/2,IJU/2,JKLOOP), &
305                     ZAPZOA2(IIU,IJU,JKLOOP)  
306   END DO
307 END IF
308 !JDJDJDJD 291196
309 ENDIF
310 !JDJDJDJD 291196
311 !
312 !-------------------------------------------------------------------------------
313 !
314 !*       4.   COMPUTE ZXHAT AND ZYHAT AT MASS POINTS
315 !              -------------------------------------
316 !
317 ZXHATM(1:IIU-1,1) = .5*(PXHAT(1:IIU-1)+PXHAT(2:IIU))
318 ZXHATM(IIU,1)     = 2.*PXHAT(IIU)-ZXHATM(IIU-1,1)
319 ZXHATM(:,2:IJU)   = SPREAD(ZXHATM(:,1),2,IJU-1)
320 !
321 ZYHATM(1,1:IJU-1) = .5*(PYHAT(1:IJU-1)+PYHAT(2:IJU))
322 ZYHATM(1,IJU)     = 2.*PYHAT(IJU)-ZYHATM(1,IJU-1)
323 ZYHATM(2:IIU,:)   = SPREAD(ZYHATM(1,:),1,IIU-1)
324 !
325 !-----------------------------------------------------------------------------
326 !
327 !*       5.   COMPUTE LATITUDES AND LONGITUDES AT MASS POINTS
328 !              -------------------------------------------------
329 !
330 CALL SM_LATLON(PLATOR,PLONOR,ZXHATM,ZYHATM,PLAT,PLON)
331 !
332 !-----------------------------------------------------------------------------
333 !
334 !*       6.  COMPUTE  MAP FACTOR AT MASS POINTS
335 !             -----------------------------------
336 !
337   IF (XRPK<0.) THEN     ! projection from north pole 
338     ZRPK=-XRPK
339     ZLAT0=-XLAT0
340     ZLAT(:,:)=-PLAT(:,:)
341   ELSE                  ! projection from south pole
342     ZRPK=XRPK
343     ZLAT0=XLAT0
344     ZLAT(:,:)=PLAT(:,:)
345   ENDIF    
346 !
347 ZRDSDG = XPI/180.
348 ZCLAT0 = COS(ZRDSDG*ZLAT0)
349 ZSLAT0 = SIN(ZRDSDG*ZLAT0)
350 !
351 IF ((ABS(ZRPK-1.)>1.E-10).AND. (ANY(ABS(COS(ZRDSDG*ZLAT))<1.E-10))) THEN
352   WRITE(ILUOUT,*) 'Error in SM_GRIDPROJ : '
353   WRITE(ILUOUT,*) 'pole in the domain, but not with stereopolar projection'
354   STOP
355 ENDIF
356 !
357 IF (ABS(ZCLAT0)<1.E-10 .AND. (ABS(ZRPK-1.)<1.E-10)) THEN
358   PMAP(:,:) = (1.+ZSLAT0)/(1.+SIN(ZRDSDG*ZLAT(:,:)))
359 ELSE
360   WHERE (ABS(COS(ZRDSDG*ZLAT(:,:)))>1.E-10)
361     PMAP(:,:) = ((ZCLAT0/COS(ZRDSDG*ZLAT(:,:)))**(1.-ZRPK))      &
362               * ((1.+ZSLAT0)/(1.+SIN(ZRDSDG*ZLAT(:,:))))**ZRPK
363   ELSEWHERE
364     PMAP(:,:) = (1.+ZSLAT0)/(1.+SIN(ZRDSDG*ZLAT(:,:)))
365   ENDWHERE
366 END IF
367 !
368 IF(NVERB >= 10) THEN                               !Value control
369   WRITE(ILUOUT,*) 'Some PMAP values:'
370   WRITE(ILUOUT,*) PMAP(1,1),PMAP(IIU/2,IJU/2),PMAP(IIU,IJU)  
371 END IF
372 !
373 !-------------------------------------------------------------------------------
374 !
375 !*       7.   COMPUTE LOCAL MESH-SIZES AT MASS POINTS
376 !              --------------------------------------
377 !
378 PDXHAT(1:IIU-1)  = PXHAT(2:IIU) - PXHAT(1:IIU-1)
379 PDXHAT(IIU)      = PDXHAT(IIU-1)
380 !
381 PDYHAT(1:IJU-1)  = PYHAT(2:IJU) - PYHAT(1:IJU-1)
382 PDYHAT(IJU)      = PDYHAT(IJU-1)
383 !
384 !JDJDJDJD 291196
385               print*,'CCONF=',CCONF,' CSTORAGE_TYPE=',CSTORAGE_TYPE
386 IF((CCONF /= 'POSTP') .OR. (CCONF == 'POSTP' .AND. CSTORAGE_TYPE /= 'PG' &
387                                              .AND. CSTORAGE_TYPE /= 'SU' ))THEN
388 !JDJDJDJD 291196
389 ZDZ(:,:,1:IKU-1) = PZZ(:,:,2:IKU) - PZZ(:,:,1:IKU-1)
390 ZDZ(:,:,IKU)     = ZDZ(:,:,IKU-1)
391 !
392 !-------------------------------------------------------------------------------
393 !
394 !*       8.    COMPUTE J AT MASS POINTS
395 !              -------------------------
396 !
397 PJ(:,:,:)  =  ZAPZOA2(:,:,:)                                                   & 
398            * SPREAD(                                                           &
399             (1/PMAP(:,:)**2)*(SPREAD(PDXHAT(:),2,IJU)*SPREAD(PDYHAT(:),1,IIU)) &
400              ,3,IKU) * ZDZ(:,:,:) 
401 !JDJDJDJD 291196
402 ENDIF
403 !JDJDJDJD 291196
404 !
405
406
407 RETURN
408 !-----------------------------------------------------------------------------
409 END SUBROUTINE SM_GRIDPROJ
410 !-----------------------------------------------------------------------------
411 !
412 !-----------------------------------------------------------------------------
413 !
414 !
415 !
416 !*              2.   ROUTINE SM_LATLON_S  (Scalar Version)
417 !                    -------------------
418 !----------------------------------------------------------------------------
419 !      #################################################
420        SUBROUTINE SM_LATLON_S(PLATOR,PLONOR,PXHATM,PYHATM,PLAT,PLON)
421 !      #################################################
422 !
423 !!****  *SM_LATLON_S * - Routine to compute geographical coordinates
424 !!
425 !!     PURPOSE
426 !!     -------
427 !        This routine computes the latitude and longitude of
428 !      a single point from  the cartesian conformal coordinates
429 !        Five map projections are available: 
430 !      - polar-stereographic from south pole  (XRPK=1),
431 !      - lambert conformal from south pole  (0<XRPK<1),
432 !      - mercator                             (XRPK=0),
433 !      - lambert conformal from north pole (-1<XRPK<0),
434 !      - polar-stereographic from north pole  (XRPK=-1).
435 !
436 !
437 !!**   METHOD
438 !!     ------
439 !!       Spherical earth approximation is used. Longitude origin is 
440 !!     set in Greenwich, and is positive eastwards. An anticlockwise 
441 !!     rotation of XBETA degrees is applied to the conformal frame 
442 !!     with respect to the geographical directions.
443 !!
444 !!       WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
445 !!
446 !!     EXTERNAL
447 !!     --------
448 !!       None
449 !!
450 !!     EXPLICIT ARGUMENTS
451 !!     ------------------
452 !!       PXHAT,PYHAT(:)  : 1D arrays of the "velocity" gridpoints
453 !!                         cartesian conformal coordinates (meters,input).
454 !!       PLATOR   : Latitude of the (1,1) point of the "mass" grid
455 !!                      (degrees,input);
456 !!       PLONOR   : Longitude of the (1,1) point of the "mass" grid
457 !!                      (degrees,input);
458 !!       PXHATM   : conformal coordinate x  (meters, mass-grid, input)
459 !!       PYHATM   : conformal coordinate y  (meters, mass-grid, input)
460 !!       PLAT     : latitude                (degrees, mass-grid, output)
461 !!       PLON     : longitude               (degrees, mass-grid, output)
462 !!
463 !!     IMPLICIT ARGUMENTS
464 !!     ------------------
465 !!       Module MODD_CST        : contains Physical constants  
466 !!          XPI        : Pi;    
467 !!          XRADIUS    : Earth radius (meters);
468 !!
469 !!       Module MODD_GRID       : contains spatial grid variables
470 !!          XLON0,XLAT0  : Reference latitude and longitude for 
471 !!                            the conformal projection (degrees);
472 !!          XBETA        : Rotation angle of the conformal frame
473 !!                            with respect to the geographical  
474 !!                            north (degrees);
475 !!          XRPK         : Projection constant (0 Mercator,
476 !!                            0<XRPK<1 Lambert, 1 Polar-stereographic)
477 !!
478 !!     REFERENCE
479 !!     ---------
480 !!      Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
481 !!            commun CNRM-LA, specifications techniques", 
482 !!            Note CNRM/GMME, 26, 139p, (Chapter 2).
483 !!      Ducrocq V., 1994, "Generation de la grille dans le modele",
484 !!            Note interne MNH, 5 mai, 3p.
485 !!      Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
486 !!            Internal note ARPEGE/ALADIN, february 27,28p.
487 !!      Levallois J., 1970, "Geodesie generale", Tome 2, Collection
488 !!             de l'IGN, Eyrolles, Paris, 408p.
489 !!       
490 !!     AUTHOR
491 !!     ------
492 !!      P.M.       *LA*
493 !!
494 !!     MODIFICATION
495 !!     ------------
496 !!       Original  PM  24/05/94
497 !!       Updated   PM  27/07/94
498 !!       Updated   VD  23/08/94
499 !!       Updated   VM  24/10/95 projection from north pole (XRPK<0) and 
500 !!                              longitudes set between XLON0-180. and XLON0+180.
501 !!
502 !-------------------------------------------------------------------------------
503 !
504 !*     0.     DECLARATIONS
505 !             ------------
506 !
507 USE MODD_CST
508 USE MODD_GRID         
509 !
510 IMPLICIT NONE
511 !
512 !*     0.1    Declarations of arguments and results
513 !
514 REAL,               INTENT(IN) :: PLATOR ! Latitude of the origine point
515 REAL,               INTENT(IN) :: PLONOR ! Longitude of the origine point
516 REAL,               INTENT(IN) :: PXHATM,PYHATM ! given conformal coordinates of the 
517                                                 ! proccessed point (meters);
518 REAL,               INTENT(OUT):: PLAT,PLON ! returned geographic latitude and 
519                                             ! longitude of the processed point 
520                                             ! (degrees).
521 !
522 !*     0.2    Declarations of local variables
523
524 REAL :: ZRPK,ZBETA,ZLAT0,ZLON0,ZLATOR,ZLONOR,ZYHATM
525 REAL :: ZRDSDG,ZCLAT0,ZSLAT0,ZCLATOR,ZSLATOR
526 REAL :: ZXBM0,ZYBM0,ZRO0,ZGA0 
527 !! JDJDJDJDJD Modif pour supporter des calculs intermediaires de capacite>32bits
528 !REAL :: ZXP,ZYP,ZEPSI,ZT1,ZCGAM,ZSGAM,ZRACLAT0
529 REAL :: ZXP,ZYP,ZEPSI,ZCGAM,ZSGAM,ZRACLAT0
530 REAL(KIND=8) :: ZT1
531 !
532 !REAL :: ZATA,ZRO2,ZT2,ZXMI0,ZYMI0
533 REAL :: ZATA,ZRO2,ZXMI0,ZYMI0,ZJD3
534 REAL(KIND=8) :: ZT2,ZJD1
535 !!!! JDJDJD
536 !
537 !--------------------------------------------------------------------------------
538 !
539 !*     1.     PRELIMINARY CALCULATIONS FOR ALL PROJECTIONS
540 !             --------------------------------------------
541 !
542 ZRDSDG = XPI/180.             ! Degree to radian conversion factor
543 ZEPSI  = 10.*EPSILON(1.)      ! A small number
544 !
545 ! By definition, (PLONOR,PLATOR) are the geographical 
546 ! coordinates, and (ZXBM0,ZYBM0) the conformal cartesian 
547 !! coordinates of the (1,1) point of the "mass" grid.
548 ! coordinates x=0, y=0 of the grid.
549 !
550 ZXBM0 = 0.
551 ZYBM0 = 0.
552
553 !
554 !--------------------------------------------------------------------------------
555 !
556 !*     2.     POLAR STEREOGRAPHIC AND LAMBERT CONFORMAL CASES
557 !             -----------------------------------------------
558 !                   (XRPK=1 P-stereo, 0<XRPK<1 Lambert)
559 !
560 IF (XRPK /= 0.) THEN
561 !
562   IF (XRPK<0.) THEN     ! projection from north pole
563     ZRPK=-XRPK
564     ZBETA=-XBETA
565     ZLAT0=-XLAT0
566     ZLON0=XLON0+180.
567     ZLATOR=-PLATOR
568     ZLONOR=PLONOR+180.
569     ZYHATM=-PYHATM
570     ZYBM0=-ZYBM0
571   ELSE                  ! projection from south pole
572     ZRPK=XRPK
573     ZBETA=XBETA
574     ZLAT0=XLAT0
575     ZLON0=XLON0
576     ZLATOR=PLATOR
577     ZLONOR=PLONOR
578     ZYHATM=PYHATM
579   ENDIF    
580 !
581 !
582 !*     2.1    Preliminary calculations
583 !
584   ZCLAT0  = COS(ZRDSDG*ZLAT0)
585   ZSLAT0  = SIN(ZRDSDG*ZLAT0)
586   ZCLATOR = COS(ZRDSDG*ZLATOR)
587   ZSLATOR = SIN(ZRDSDG*ZLATOR)
588   ZRO0    = (XRADIUS/ZRPK)*(ABS(ZCLAT0))**(1.-ZRPK)     &
589           * ((1.+ZSLAT0)*ABS(ZCLATOR)/(1.+ZSLATOR))**ZRPK
590   ZGA0    = (ZRPK*(ZLONOR-ZLON0)-ZBETA)*ZRDSDG
591   ZXP     = ZXBM0-ZRO0*SIN(ZGA0)
592   ZYP     = ZYBM0+ZRO0*COS(ZGA0)
593 !
594 !*    2.2    Longitude
595 !
596   IF(ABS(ZYHATM-ZYP) < ZEPSI.AND.ABS(PXHATM-ZXP) < ZEPSI)THEN
597     ZATA = 0.
598   ELSE
599     ZATA = ATAN2(-(ZXP-PXHATM),(ZYP-ZYHATM))/ZRDSDG
600   ENDIF
601   !
602   PLON = (ZBETA+ZATA)/ZRPK+ZLON0
603 !
604 !*   2.3     Latitude
605 !
606   ZRO2 = (PXHATM-ZXP)**2+(ZYHATM-ZYP)**2
607 !! JDJDJDJDJD Modif pour supporter des calculs intermediaires de capacite>32bits
608   ZJD1 = XRADIUS*(ABS(ZCLAT0))**(1.-ZRPK)
609   ZT1  = (ZJD1)**(2./ZRPK)   &
610        * (1+ZSLAT0)**2
611 ! ZT1  = (XRADIUS*(ABS(ZCLAT0))**(1.-ZRPK))**(2./ZRPK)   &
612 !      * (1+ZSLAT0)**2
613   ZJD3 = (ZRPK**2*ZRO2)
614   ZT2  = ZJD3
615   ZT2 = ZT2**(1./ZRPK)
616 ! ZT2  = (ZRPK**2*ZRO2)**(1./ZRPK)
617   !
618   ZJD1 = (ZT1-ZT2)/(ZT1+ZT2)
619   ZJD1 = ACOS(ZJD1)
620   ZJD3 = ZJD1
621   PLAT = (XPI/2.-ZJD3)/ZRDSDG
622 ! PLAT = (XPI/2.-ACOS((ZT1-ZT2)/(ZT1+ZT2)))/ZRDSDG
623 !! JDJDJDJDJD
624 !
625   IF (XRPK<0.) THEN     ! projection from north pole 
626     PLAT=-PLAT
627     PLON=PLON-180.
628   ENDIF    
629 !
630 !---------------------------------------------------------------------------------
631 !
632 !*  3.        MERCATOR PROJECTION WITH ROTATION
633 !             ---------------------------------
634 !                       (XRPK=0)
635 !
636 ELSE
637 !
638 !*  3.1       Preliminary calculations
639 !
640   ZCGAM    = COS(-ZRDSDG*XBETA)
641   ZSGAM    = SIN(-ZRDSDG*XBETA)
642   ZRACLAT0 = XRADIUS*COS(ZRDSDG*XLAT0)
643 !
644 !*  3.2       Longitude
645 !
646   ZXMI0 = PXHATM-ZXBM0
647   ZYMI0 = PYHATM-ZYBM0
648   !
649   PLON  = (ZXMI0*ZCGAM+ZYMI0*ZSGAM)/(ZRACLAT0*ZRDSDG)+PLONOR
650 !
651 !*  3.3       Latitude
652 !
653   ZT1  = LOG(TAN(XPI/4.+PLATOR*ZRDSDG/2.))
654   ZT2  = (-ZXMI0*ZSGAM+ZYMI0*ZCGAM)/ZRACLAT0
655   !
656   PLAT = (-XPI/2.+2.*ATAN(EXP(ZT1+ZT2)))/ZRDSDG
657 !
658 !---------------------------------------------------------------------------------
659 !
660 !*  4.        EXIT
661 !             ----
662 !
663 END IF
664 PLON=PLON+NINT((XLON0-PLON)/360.)*360.
665 RETURN
666 !--------------------------------------------------------------------------------
667 END SUBROUTINE SM_LATLON_S
668 !-------------------------------------------------------------------------------
669 !
670 !---------------------------------------------------------------------------------
671 !
672 !*              3.   ROUTINE SM_LATLON_A  (Array Version )
673 !                    -------------------
674 !--------------------------------------------------------------------------------
675 !      ###################################################
676        SUBROUTINE SM_LATLON_A(PLATOR,PLONOR,  &
677                               PXHATM,PYHATM,PLAT,PLON)
678 !      ###################################################
679 !
680 !!****  *SM_LATLON_A * - Routine to compute geographical coordinates
681 !!
682 !!     PURPOSE
683 !!     -------
684 !        This routine computes the latitude and longitude of
685 !      an array given in cartesian conformal coordinates
686 !        Five map projections are available: 
687 !      - polar-stereographic from south pole  (XRPK=1),
688 !      - lambert conformal from south pole  (0<XRPK<1),
689 !      - mercator                             (XRPK=0),
690 !      - lambert conformal from north pole (-1<XRPK<0),
691 !      - polar-stereographic from north pole  (XRPK=-1).
692 !
693 !
694 !!**   METHOD
695 !!     ------
696 !!       Spherical earth approximation is used. Longitude origin is 
697 !!     set in Greenwich, and is positive eastwards. An anticlockwise 
698 !!     rotation of XBETA degrees is applied to the conformal frame 
699 !!     with respect to the geographical directions.
700 !!
701 !!       WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
702 !!
703 !!     EXTERNAL
704 !!     --------
705 !!       None
706 !!
707 !!     EXPLICIT ARGUMENTS
708 !!     ------------------
709 !!       PXHAT,PYHAT(:)  : 1D arrays of the "velocity" gridpoints
710 !!                         cartesian conformal coordinates (meters,input).
711 !!       PLATOR   : Latitude of the (1,1) point of the "mass" grid
712 !!                      (degrees,input);
713 !!       PLONOR   : Longitude of the (1,1) point of the "mass" grid
714 !!                      (degrees,input);
715 !!       PXHATM   : conformal coordinate x  (meters, mass-grid, input)
716 !!       PYHATM   : conformal coordinate y  (meters, mass-grid, input)
717 !!       PLAT    : latitude                (degrees, mass-grid, output)
718 !!       PLON    : longitude               (degrees, mass-grid, output)
719 !!
720 !!
721 !!     IMPLICIT ARGUMENTS
722 !!     ------------------
723 !!       Module MODD_CST      : contains Physical constants
724 !!          XPI           : Pi;    
725 !!          XRADIUS       : Earth radius (meters);
726 !!
727 !!       Module MODD_GRID     : contains spatial grid variables
728 !!          XLON0,XLAT0   : Reference latitude and longitude for 
729 !!                          the conformal projection (degrees);
730 !!          XBETA         : Rotation angle of the conformal frame
731 !!                          with respect to the geographical  
732 !!                          north (degrees);
733 !!          XRPK          : Projection constant (0 Mercator,
734 !!                          0<XRPK<1 Lambert, 1 Polar-stereographic);
735 !!
736 !!     REFERENCE
737 !!     ---------
738 !!      Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
739 !!            commun CNRM-LA, specifications techniques", 
740 !!            Note CNRM/GMME, 26, 139p, (Chapter 2).
741 !!      Ducrocq V., 1994, "Generation de la grille dans le modele",
742 !!            Note interne MNH, 5 mai, 3p.
743 !!      Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
744 !!            Internal note ARPEGE/ALADIN, february 27,28p.
745 !!      Levallois J., 1970, "Geodesie generale", Tome 2, Collection
746 !!             de l'IGN, Eyrolles, Paris, 408p.
747 !!       
748 !!     AUTHOR
749 !!     ------
750 !!      P.M.       *LA*
751 !!
752 !!     MODIFICATION
753 !!     ------------
754 !!       Original  PM  24/05/94
755 !!       Updated   PM  27/07/94
756 !!       Updated   VD  23/08/94
757 !!       Updated   VM  24/10/95 projection from north pole (XRPK<0) and 
758 !!                              longitudes set between XLON0-180. and XLON0+180.
759 !!
760 !-------------------------------------------------------------------------------
761 !
762 !*     0.     DECLARATIONS
763 !             ------------
764 !
765 USE MODD_CST
766 USE MODD_GRID              
767 !
768 IMPLICIT NONE
769 !
770 !*     0.1    Declarations of arguments and results
771 !
772 REAL,                 INTENT(IN) :: PLATOR ! Latitude of the origine point
773 REAL,                 INTENT(IN) :: PLONOR ! Longitude of the origine point
774 REAL, DIMENSION(:,:), INTENT(IN) :: PXHATM,PYHATM   
775                                 ! given conformal coordinates of the 
776                                 ! processed points (meters);
777 REAL, DIMENSION(:,:), INTENT(OUT):: PLAT,PLON    
778                                 ! returned geographic latitudes and 
779                                 ! longitudes of the processed points 
780                                 ! (degrees).
781 !
782 !*     0.2    Declarations of local variables
783
784 REAL, DIMENSION(SIZE(PYHATM,1),SIZE(PYHATM,2)) :: ZYHATM
785 REAL :: ZRPK,ZBETA,ZLAT0,ZLON0,ZLATOR,ZLONOR
786 REAL :: ZRDSDG,ZCLAT0,ZSLAT0,ZCLATOR,ZSLATOR
787 REAL :: ZXBM0,ZYBM0,ZRO0,ZGA0 
788 !! JDJDJDJDJD Modif pour supporter des calculs intermediaires de capacite>32bits
789 !REAL :: ZXP,ZYP,ZEPSI,ZT1,ZCGAM,ZSGAM,ZRACLAT0
790 REAL :: ZXP,ZYP,ZEPSI,ZCGAM,ZSGAM,ZRACLAT0
791 REAL(KIND=8) :: ZT1,ZJD4,ZJD5
792 REAL :: ZRPK2
793 !!! JDJDJDJDJD 
794 !
795 !! JDJDJDJDJD Modif pour supporter des calculs intermediaires de capacite>32bits
796 !REAL, DIMENSION(SIZE(PXHATM,1),SIZE(PXHATM,2)) :: ZATA,ZRO2,ZT2,ZXMI0,ZYMI0
797 REAL, DIMENSION(SIZE(PXHATM,1),SIZE(PXHATM,2)) :: ZATA,ZRO2,ZXMI0,ZYMI0,ZJD3
798 REAL(KIND=8), DIMENSION(SIZE(PXHATM,1),SIZE(PXHATM,2)) :: ZT2,ZJD1,ZJD2
799 !!! JDJDJDJDJD 
800 !
801 !--------------------------------------------------------------------------------
802 !
803 !*     1.     Preliminary calculations for all projections
804 !             --------------------------------------------
805 !
806 ZRDSDG = XPI/180.         ! Degree to radian conversion factor
807 ZEPSI  = 10.*EPSILON(1.)      ! A small number
808 !
809 ! By definition, (PLONOR,PLATOR) are the geographical 
810 ! coordinates, and (ZXBM0,ZYBM0) the conformal cartesian 
811 ! coordinates x=0, y=0 of the grid.
812 !! coordinates of the (1,1) point of the "mass" grid.
813 !
814 ZXBM0 = 0.
815 ZYBM0 = 0.
816 !
817 !-------------------------------------------------------------------------------
818 !
819 !*     2.     POLAR STEREOGRAPHIC AND LAMBERT CONFORMAL CASES
820 !             -----------------------------------------------
821 !                   (XRPK=1 P-stereo, 0<XRPK<1 Lambert)
822 !
823 IF(XRPK /= 0.) THEN
824 !
825   IF (XRPK<0.) THEN     ! projection from north pole
826     ZRPK=-XRPK
827     ZBETA=-XBETA
828     ZLAT0=-XLAT0
829     ZLON0=XLON0+180.
830     ZLATOR=-PLATOR
831     ZLONOR=PLONOR+180.
832     ZYHATM(:,:)=-PYHATM(:,:)
833     ZYBM0=-ZYBM0
834   ELSE                  ! projection from south pole
835     ZRPK=XRPK
836     ZBETA=XBETA
837     ZLAT0=XLAT0
838     ZLON0=XLON0
839     ZLATOR=PLATOR
840     ZLONOR=PLONOR
841     ZYHATM(:,:)=PYHATM(:,:)
842   ENDIF    
843 !
844 !*     2.1    Preliminary calculations
845 !
846   ZCLAT0  = COS(ZRDSDG*ZLAT0)
847   ZSLAT0  = SIN(ZRDSDG*ZLAT0)
848   ZCLATOR = COS(ZRDSDG*ZLATOR)
849   ZSLATOR = SIN(ZRDSDG*ZLATOR)
850   ZRO0    = (XRADIUS/ZRPK)*(ABS(ZCLAT0))**(1.-ZRPK)     &
851           * ((1.+ZSLAT0)*ABS(ZCLATOR)/(1.+ZSLATOR))**ZRPK
852   ZGA0    = (ZRPK*(ZLONOR-ZLON0)-ZBETA)*ZRDSDG
853   ZXP     = ZXBM0-ZRO0*SIN(ZGA0)
854   ZYP     = ZYBM0+ZRO0*COS(ZGA0)
855 !
856 !*    2.2    Longitude
857 !
858   WHERE  (ABS(ZYHATM(:,:)-ZYP) < ZEPSI    &
859      .AND.ABS(PXHATM(:,:)-ZXP) < ZEPSI)
860     ZATA(:,:) = 0.
861   ELSEWHERE
862     ZATA(:,:) = ATAN2(-(ZXP-PXHATM(:,:)),(ZYP-ZYHATM(:,:)))/ZRDSDG
863   END WHERE
864   !
865   PLON(:,:) = (ZBETA+ZATA(:,:))/ZRPK+ZLON0
866 !
867 !*   2.3     Latitude
868 !
869   ZRO2(:,:) = (PXHATM(:,:)-ZXP)**2+(ZYHATM(:,:)-ZYP)**2
870   ZJD4      = (XRADIUS*(ABS(ZCLAT0))**(1.-ZRPK))
871   ZJD5      = ZJD4**(2./ZRPK)
872   ZT1       = ZJD5 * (1+ZSLAT0)**2
873 ! ZT1       = (XRADIUS*(ABS(ZCLAT0))**(1.-ZRPK))**(2./ZRPK)   &
874 !          * (1+ZSLAT0)**2
875   ZRPK2 = ZRPK**2
876   ZJD3(:,:) = (ZRPK2*ZRO2(:,:))
877   ZT2(:,:)  = ZJD3(:,:)
878   ZT2(:,:)  = ZT2(:,:)**(1./ZRPK)
879 ! ZT2(:,:)  = (ZRPK**2*ZRO2(:,:))**(1./ZRPK)
880 !
881 !! JDJDJDJDJD Modif pour supporter des calculs intermediaires de capacite>32bits
882   ZJD1(:,:) = (ZT1-ZT2(:,:))
883   ZJD2(:,:) = (ZT1+ZT2(:,:))
884   ZJD1(:,:) = ZJD1(:,:)/ZJD2(:,:)
885   ZJD1(:,:) = ACOS(ZJD1(:,:))
886   ZJD3(:,:) = ZJD1(:,:)
887   PLAT(:,:) = (XPI/2.-ZJD3(:,:))/ZRDSDG
888 ! PLAT(:,:) = (XPI/2.-ACOS((ZT1-ZT2(:,:))/(ZT1+ZT2(:,:))))/ZRDSDG
889 !! JDJDJDJDJD 
890 !
891   IF (XRPK<0.) THEN     ! projection from north pole
892     PLAT(:,:)=-PLAT(:,:)
893     PLON(:,:)=PLON(:,:)+180.
894   ENDIF
895 !
896 !-------------------------------------------------------------------------------
897 !
898 !*  3.        MERCATOR PROJECTION WITH ROTATION
899 !             ---------------------------------
900 !                       (XRPK=0)
901 !
902 ELSE
903 !
904 !*  3.1       Preliminary calculations
905 !
906   ZCGAM    = COS(-ZRDSDG*XBETA)
907   ZSGAM    = SIN(-ZRDSDG*XBETA)
908   ZRACLAT0 = XRADIUS*COS(ZRDSDG*XLAT0)
909 !
910 !*  3.2       Longitude
911 !
912   ZXMI0(:,:) = PXHATM(:,:)-ZXBM0
913   ZYMI0(:,:) = PYHATM(:,:)-ZYBM0
914   !
915   PLON(:,:) = (ZXMI0(:,:)*ZCGAM+ZYMI0(:,:)*ZSGAM)     &
916             / (ZRACLAT0*ZRDSDG)+PLONOR
917 !
918 !*  3.3       Latitude
919 !
920   ZT1       = ALOG(TAN(XPI/4.+PLATOR*ZRDSDG/2.))
921   ZT2(:,:)  = (-ZXMI0(:,:)*ZSGAM+ZYMI0(:,:)*ZCGAM)/ZRACLAT0
922   !
923   PLAT(:,:) = (-XPI/2.+2.*ATAN(EXP(ZT1+ZT2(:,:))))/ZRDSDG
924 !
925 !---------------------------------------------------------------------------------
926 !
927 !*  4.        EXIT
928 !             ----
929 !
930 END IF
931 PLON(:,:)=PLON(:,:)+NINT((XLON0-PLON(:,:))/360.)*360.
932 RETURN
933 !---------------------------------------------------------------------------------
934 END SUBROUTINE SM_LATLON_A
935 !---------------------------------------------------------------------------------
936 !
937 !---------------------------------------------------------------------------------
938 !
939 !*              4.   ROUTINE SM_XYHAT_S (Scalar Version )
940 !                    ------------------
941 !--------------------------------------------------------------------------------
942 !      ##################################################
943        SUBROUTINE SM_XYHAT_S(PLATOR,PLONOR,  &
944                              PLAT,PLON,PXHATM,PYHATM)
945 !      ##################################################
946 !
947 !!****  *SM_XYHAT_S * - Routine to compute conformal coordinates
948 !!
949 !!     PURPOSE
950 !!     -------
951 !        This routine computes the cartesian conformal coordinates 
952 !      of a single point from  its latitude and longitude
953 !        Five map projections are available: 
954 !      - polar-stereographic from south pole  (XRPK=1),
955 !      - lambert conformal from south pole  (0<XRPK<1),
956 !      - mercator                             (XRPK=0),
957 !      - lambert conformal from north pole (-1<XRPK<0),
958 !      - polar-stereographic from north pole  (XRPK=-1).
959 !
960 !
961 !!**   METHOD
962 !!     ------
963 !!       Spherical earth approximation is used. Longitude origin is 
964 !!     set in Greenwich, and is positive eastwards. An anticlockwise 
965 !!     rotation of XBETA degrees is applied to the conformal frame 
966 !!     with respect to the geographical directions.
967 !!
968 !!       WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
969 !!
970 !!     EXTERNAL
971 !!     --------
972 !!       None
973 !!
974 !!     EXPLICIT ARGUMENTS
975 !!     ------------------
976 !!       PLATOR   : Latitude of the (1,1) point of the "mass" grid
977 !!                      (degrees,input);
978 !!       PLONOR   : Longitude of the (1,1) point of the "mass" grid
979 !!                      (degrees,input);
980 !!       PXHATM   : conformal coordinate x  (meters, mass-grid, input)
981 !!       PYHATM   : conformal coordinate y  (meters, mass-grid, input)
982 !!       PLAT    : latitude                (degrees, mass-grid, output)
983 !!       PLON    : longitude               (degrees, mass-grid, output)
984 !!
985 !!     IMPLICIT ARGUMENTS
986 !!     ------------------
987 !!       Module MODD_CST     : contains Physical constants
988 !!          XPI          : Pi;    
989 !!          XRADIUS      : Earth radius (meters);
990 !!
991 !!       Module MODD_GRID    : contains spatial grid variables
992 !!          XLON0,XLAT0  : Reference latitude and longitude for 
993 !!                         the conformal projection (degrees);
994 !!          XBETA        : Rotation angle of the conformal frame
995 !!                         with respect to the geographical  
996 !!                         north (degrees);
997 !!          XRPK         : Projection constant (0 Mercator,
998 !!                         0<XRPK<1 Lambert, 1 Polar-stereographic);
999 !!
1000 !!     REFERENCE
1001 !!     ---------
1002 !!      Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
1003 !!            commun CNRM-LA, specifications techniques", 
1004 !!            Note CNRM/GMME, 26, 139p, (Chapter 2).
1005 !!      Ducrocq V., 1994, "Generation de la grille dans le modele",
1006 !!            Note interne MNH, 5 mai, 3p.
1007 !!      Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
1008 !!            Internal note ARPEGE/ALADIN, february 27,28p.
1009 !!      Levallois J., 1970, "Geodesie generale", Tome 2, Collection
1010 !!             de l'IGN, Eyrolles, Paris, 408p.
1011 !!       
1012 !!     AUTHOR
1013 !!     ------
1014 !!      P.M.       *LA*
1015 !!
1016 !!     MODIFICATION
1017 !!     ------------
1018 !!       Original  PM  24/05/94
1019 !!       Updated   PM  27/07/94
1020 !!       Updated   VD  23/08/94
1021 !!       Updated   VM  24/10/95 projection from north pole (XRPK<0) and 
1022 !!                              longitudes set between XLON0-180. and XLON0+180.
1023 !!
1024 !-------------------------------------------------------------------------------
1025 !
1026 !*     0.     DECLARATIONS
1027 !             ------------
1028 !
1029 USE MODD_CST
1030 USE MODD_GRID              
1031 !
1032 IMPLICIT NONE
1033 !
1034 !*     0.1    Declarations of arguments and results
1035 !
1036 REAL,               INTENT(IN) :: PLATOR ! Latitude of the origine point
1037 REAL,               INTENT(IN) :: PLONOR ! Longitude of the origine point
1038 REAL,               INTENT(IN) :: PLAT,PLON 
1039                                          ! given geographic latitude and 
1040                                          ! longitude of the processed point 
1041                                          ! (degrees).
1042 REAL,               INTENT(OUT):: PXHATM,PYHATM 
1043                                          ! returned conformal coordinates of 
1044                                          ! the processed point (meters);
1045 !
1046 !*     0.2    Declarations of local variables
1047
1048 REAL :: ZRPK,ZBETA,ZLAT0,ZLON0,ZLATOR,ZLONOR
1049 REAL :: ZLAT,ZLON
1050 REAL :: ZRDSDG,ZCLAT0,ZSLAT0,ZCLATOR,ZSLATOR
1051 REAL :: ZXBM0,ZYBM0,ZRO0,ZGA0 
1052 REAL :: ZXP,ZYP,ZCGAM,ZSGAM,ZRACLAT0,ZXE,ZYE
1053 !
1054 REAL :: ZCLAT,ZSLAT,ZRO,ZGA,ZXPR,ZYPR
1055 !
1056 !--------------------------------------------------------------------------------
1057 !
1058 !*     1.     PRELIMINARY CALCULATION FOR ALL PROJECTIONS
1059 !             -------------------------------------------
1060 !
1061 ZRDSDG = XPI/180.         ! Degree to radian conversion factor
1062 !
1063 ! By definition, (PLONOR,PLATOR) are the geographical 
1064 ! coordinates of the x=0, y=0 point.
1065 !
1066 ZXBM0 = 0.
1067 ZYBM0 = 0.
1068 !
1069 ZLON=PLON
1070 ZLON=ZLON+NINT((XLON0-ZLON)/360.)*360.
1071 !
1072 ZLONOR=PLONOR
1073 ZLONOR=ZLONOR+NINT((XLON0-ZLONOR)/360.)*360.
1074 !---------------------------------------------------------------------------------
1075 !
1076 !*     2.     POLAR STEREOGRAPHIC AND LAMBERT CONFORMAL CASES
1077 !             -----------------------------------------------
1078 !                   (XRPK=1 P-stereo, 0<XRPK<1 Lambert)
1079 !
1080 IF(XRPK /= 0.) THEN
1081 !
1082   IF (XRPK<0.) THEN     ! projection from north pole
1083     ZRPK=-XRPK
1084     ZBETA=-XBETA
1085     ZLAT0=-XLAT0
1086     ZLON0=XLON0+180.
1087     ZLATOR=-PLATOR
1088     ZLONOR=ZLONOR+180.
1089     ZLAT=-PLAT
1090     ZLON=ZLON+180.
1091     ZYBM0=-ZYBM0
1092   ELSE                  ! projection from south pole
1093     ZRPK=XRPK
1094     ZBETA=XBETA
1095     ZLAT0=XLAT0
1096     ZLON0=XLON0
1097     ZLATOR=PLATOR
1098     ZLONOR=ZLONOR
1099     ZLAT=PLAT
1100     ZLON=ZLON
1101   ENDIF    
1102 !
1103 !*     2.1    Preliminary calculations
1104 !
1105   ZCLAT0  = COS(ZRDSDG*ZLAT0)
1106   ZSLAT0  = SIN(ZRDSDG*ZLAT0)
1107   ZCLATOR = COS(ZRDSDG*ZLATOR)
1108   ZSLATOR = SIN(ZRDSDG*ZLATOR)
1109   ZRO0    = (XRADIUS/ZRPK)*(ABS(ZCLAT0))**(1.-ZRPK)     &
1110           * ((1.+ZSLAT0)*ABS(ZCLATOR)/(1.+ZSLATOR))**ZRPK
1111   ZGA0    = (ZRPK*(ZLONOR-ZLON0)-ZBETA)*ZRDSDG
1112   ZXP     = ZXBM0-ZRO0*SIN(ZGA0)
1113   ZYP     = ZYBM0+ZRO0*COS(ZGA0)
1114 !
1115 !*    2.2    Conformal coordinates in meters
1116 !
1117   ZCLAT  = COS(ZRDSDG*ZLAT)
1118   ZSLAT  = SIN(ZRDSDG*ZLAT)
1119   ZRO    = (XRADIUS/ZRPK)*(ABS(ZCLAT0))**(1.-ZRPK)    &
1120          * ((1.+ZSLAT0)*ABS(ZCLAT)/(1.+ZSLAT))**ZRPK
1121   ZGA    = (ZRPK*(ZLON-ZLON0)-ZBETA)*ZRDSDG
1122 !
1123   PXHATM = ZXP+ZRO*SIN(ZGA)
1124   PYHATM = ZYP-ZRO*COS(ZGA)
1125 !
1126   IF (XRPK<0.) THEN     ! projection from north pole
1127     PYHATM=-PYHATM
1128   ENDIF
1129 !
1130 !
1131 !------------------------------------------------------------------------------
1132 !
1133 !*  3.        MERCATOR PROJECTION WITH ROTATION
1134 !             ---------------------------------
1135 !                       (XRPK=0)
1136 !
1137 ELSE
1138 !
1139 !*  3.1       Preliminary calculations
1140 !
1141   ZCGAM    = COS(-ZRDSDG*XBETA)
1142   ZSGAM    = SIN(-ZRDSDG*XBETA)
1143   ZRACLAT0 = XRADIUS*COS(ZRDSDG*XLAT0)
1144   ZXE      = ZXBM0*ZCGAM+ZYBM0*ZSGAM            &
1145            - ZRACLAT0*(PLONOR-XLON0)*ZRDSDG  
1146   ZYE      =-ZXBM0*ZSGAM+ZYBM0*ZCGAM            &
1147            - ZRACLAT0*LOG(TAN(XPI/4.+PLATOR*ZRDSDG/2.))
1148 !
1149 !*  3.2       Conformal coordinates
1150 !
1151   ZXPR   = ZRACLAT0*(ZLON-XLON0)*ZRDSDG+ZXE
1152   ZYPR   = ZRACLAT0*LOG(TAN(XPI/4.+PLAT*ZRDSDG/2.))+ZYE
1153   !
1154   PXHATM = ZXPR*ZCGAM-ZYPR*ZSGAM
1155   PYHATM = ZXPR*ZSGAM+ZYPR*ZCGAM
1156 !
1157 !-------------------------------------------------------------------------------
1158 !
1159 !*  4.        EXIT
1160 !             ----
1161 !
1162 END IF
1163 RETURN
1164 !-------------------------------------------------------------------------------
1165 END SUBROUTINE SM_XYHAT_S
1166 !-------------------------------------------------------------------------------
1167 !
1168 !-------------------------------------------------------------------------------
1169 !
1170 !*              5.   ROUTINE SM_XYHAT_A (Array Version )
1171 !                    ------------------
1172 !-------------------------------------------------------------------------------
1173 !      ################################################
1174        SUBROUTINE SM_XYHAT_A(PLATOR,PLONOR,  &
1175                              PLAT,PLON,PXHATM,PYHATM)
1176 !      ################################################
1177 !
1178 !!****  *SM_XYHAT_A * - Routine to compute conformal coordinates
1179 !!
1180 !!
1181 !!     PURPOSE
1182 !!     -------
1183 !        This routine computes the cartesian conformal coordinates 
1184 !      of an array given in latitude-longitude coordinates
1185 !        Three map projections are available: 
1186 !      - polar-stereographic (XRPK=1),
1187 !      - lambert conformal  (0<XRPK<1),
1188 !      - mercator (XRPK=0).
1189 !
1190 !
1191 !!**   METHOD
1192 !!     ------
1193 !!       Spherical earth approximation is used. Longitude origin is 
1194 !!     set in Greenwich, and is positive eastwards. An anticlockwise 
1195 !!     rotation of XBETA degrees is applied to the conformal frame 
1196 !!     with respect to the geographical directions.
1197 !!
1198 !!       WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
1199 !!
1200 !!     EXTERNAL
1201 !!     --------
1202 !!       None
1203 !!
1204 !!     EXPLICIT ARGUMENTS
1205 !!     ------------------
1206 !!       PLATOR   : Latitude of the (1,1) point of the "mass" grid
1207 !!                      (degrees,input);
1208 !!       PLONOR   : Longitude of the (1,1) point of the "mass" grid
1209 !!                      (degrees,input);
1210 !!       PXHATM   : conformal coordinate x  (meters, mass-grid, input)
1211 !!       PYHATM   : conformal coordinate y  (meters, mass-grid, input)
1212 !!       PLAT    : latitude                (degrees, mass-grid, output)
1213 !!       PLON    : longitude               (degrees, mass-grid, output)
1214 !!
1215 !!     IMPLICIT ARGUMENTS
1216 !!     ------------------
1217 !!       Module MODD_CST         : contains Physical constants
1218 !!          XPI         : Pi;    
1219 !!          XRADIUS     : Earth radius (meters);
1220 !!
1221 !!       Module MODD_GRID        : contains spatial grid variables
1222 !!          XLON0,XLAT0 : Reference latitude and longitude for 
1223 !!                        the conformal projection (degrees);
1224 !!          XBETA       : Rotation angle of the conformal frame
1225 !!                        with respect to the geographical  
1226 !!                        north (degrees);
1227 !!          XRPK        : Projection constant (0 Mercator,
1228 !!                        0<XRPK<1 Lambert, 1 Polar-stereographic);
1229 !!
1230 !!     REFERENCE
1231 !!     ---------
1232 !!      Asencio N. et al., 1994, "Le projet de modele non-hydrostatique
1233 !!            commun CNRM-LA, specifications techniques", 
1234 !!            Note CNRM/GMME, 26, 139p, (Chapter 2).
1235 !!      Ducrocq V., 1994, "Generation de la grille dans le modele",
1236 !!            Note interne MNH, 5 mai, 3p.
1237 !!      Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
1238 !!            Internal note ARPEGE/ALADIN, february 27,28p.
1239 !!      Levallois J., 1970, "Geodesie generale", Tome 2, Collection
1240 !!             de l'IGN, Eyrolles, Paris, 408p.
1241 !!       
1242 !!     AUTHOR
1243 !!     ------
1244 !!      P.M.       *LA*
1245 !!
1246 !!     MODIFICATION
1247 !!     ------------
1248 !!       Original PM  24/05/94
1249 !!       Updated  PM  27/07/94
1250 !!       Updated  VD  23/08/94
1251 !!       Updated  VM  24/10/95 projection from north pole (XRPK<0) and 
1252 !!                             longitudes set between XLON0-180. and XLON0+180.
1253 !!
1254 !-------------------------------------------------------------------------------
1255 !
1256 !*     0.     DECLARATIONS
1257 !             ------------
1258 !
1259 USE MODD_CST
1260 USE MODD_GRID       
1261 !
1262 IMPLICIT NONE
1263 !
1264 !*     0.1    Declarations of arguments and results
1265 !
1266 REAL,                INTENT(IN) :: PLATOR ! Latitude of the origine point
1267 REAL,                INTENT(IN) :: PLONOR ! Longitude of the origine point
1268 REAL, DIMENSION(:,:), INTENT(IN):: PLAT,PLON     
1269                                           ! given geographic latitude and 
1270                                           ! longitude of the processed  
1271                                           ! array (degrees).
1272 REAL, DIMENSION(:,:), INTENT(OUT):: PXHATM,PYHATM  
1273                                           ! returned conformal coordinates of 
1274                                           ! the processed array (meters);
1275 !
1276 !*     0.2    Declarations of local variables
1277
1278 REAL,DIMENSION(SIZE(PLAT,1),SIZE(PLAT,2)) :: ZLAT,ZLON
1279 REAL :: ZRPK,ZBETA,ZLAT0,ZLON0,ZLATOR,ZLONOR
1280 REAL :: ZRDSDG,ZCLAT0,ZSLAT0,ZCLATOR,ZSLATOR
1281 REAL :: ZXBM0,ZYBM0,ZRO0,ZGA0 
1282 REAL :: ZXP,ZYP,ZCGAM,ZSGAM,ZRACLAT0,ZXE,ZYE
1283 !
1284 REAL,DIMENSION(SIZE(PLAT,1),SIZE(PLAT,2)) :: ZCLAT,ZSLAT,ZRO,ZGA,ZXPR,ZYPR
1285 !
1286 !
1287 !-------------------------------------------------------------------------------
1288 !
1289 !*     1.     PRELIMINARY CALCULATION FOR ALL PROJECTIONS
1290 !             -------------------------------------------
1291 !
1292 ZRDSDG = XPI/180.         ! Degree to radian conversion factor
1293 !
1294 ! By definition, (PLONOR,PLATOR) are the geographical 
1295 ! coordinates of the x=0, y=0 point.
1296 !
1297 ZXBM0 = 0.
1298 ZYBM0 = 0.
1299 !
1300 ZLON(:,:)=PLON(:,:)
1301 ZLON(:,:)=ZLON(:,:)+NINT((XLON0-ZLON(:,:))/360.)*360.
1302 !
1303 ZLONOR=PLONOR
1304 ZLONOR=ZLONOR+NINT((XLON0-ZLONOR)/360.)*360.
1305 !------------------------------------------------------------------------------
1306 !
1307 !*     2.     POLAR SEREOGRAPHIC AND LAMBERT CONFORMAL CASES
1308 !             ----------------------------------------------
1309 !                   (XRPK=1 P-stereo, 0<XRPK<1 Lambert)
1310 !
1311 IF(XRPK  /=  0.) THEN
1312 !
1313   IF (XRPK<0.) THEN     ! projection from north pole
1314     ZRPK=-XRPK
1315     ZBETA=-XBETA
1316     ZLAT0=-XLAT0
1317     ZLON0=XLON0+180.
1318     ZLATOR=-PLATOR
1319     ZLONOR=ZLONOR+180.
1320     ZLAT(:,:)=-PLAT(:,:)
1321     ZLON(:,:)=ZLON(:,:)+180.
1322     ZYBM0=-ZYBM0
1323   ELSE                  ! projection from south pole
1324     ZRPK=XRPK
1325     ZBETA=XBETA
1326     ZLAT0=XLAT0
1327     ZLON0=XLON0
1328     ZLATOR=PLATOR
1329     ZLONOR=ZLONOR
1330     ZLAT(:,:)=PLAT(:,:)
1331     ZLON(:,:)=ZLON(:,:)
1332   ENDIF    
1333 !
1334 !*     2.1    Preliminary calculations
1335 !
1336   ZCLAT0  = COS(ZRDSDG*ZLAT0)
1337   ZSLAT0  = SIN(ZRDSDG*ZLAT0)
1338   ZCLATOR = COS(ZRDSDG*ZLATOR)
1339   ZSLATOR = SIN(ZRDSDG*ZLATOR)
1340   ZRO0    = (XRADIUS/ZRPK)*(ABS(ZCLAT0))**(1.-ZRPK)    &
1341           * ((1.+ZSLAT0)*ABS(ZCLATOR)/(1.+ZSLATOR))**ZRPK
1342   ZGA0    = (ZRPK*(ZLONOR-ZLON0)-ZBETA)*ZRDSDG
1343   ZXP     = ZXBM0-ZRO0*SIN(ZGA0)
1344   ZYP     = ZYBM0+ZRO0*COS(ZGA0)
1345 !
1346 !*    2.2    Conformal coordinates in meters
1347 !
1348   ZCLAT(:,:)  = COS(ZRDSDG*ZLAT(:,:))
1349   ZSLAT(:,:)  = SIN(ZRDSDG*ZLAT(:,:))
1350   ZRO(:,:)    = (XRADIUS/ZRPK)*(ABS(ZCLAT0))**(1.-ZRPK)    &
1351               * ((1.+ZSLAT0)*ABS(ZCLAT(:,:))/(1.+ZSLAT(:,:)))**ZRPK
1352   ZGA(:,:)    = (ZRPK*(ZLON(:,:)-ZLON0)-ZBETA)*ZRDSDG
1353 !
1354   PXHATM(:,:) = ZXP+ZRO(:,:)*SIN(ZGA(:,:))
1355   PYHATM(:,:) = ZYP-ZRO(:,:)*COS(ZGA(:,:))
1356 !
1357   IF (XRPK<0.) THEN     ! projection from north pole
1358     PYHATM(:,:)=-PYHATM(:,:)
1359   ENDIF
1360 !
1361 !-------------------------------------------------------------------------------
1362 !
1363 !*  3.        MERCATOR PROJECTION WITH ROTATION
1364 !             ---------------------------------
1365 !                       (XRPK=0)
1366 !
1367 ELSE
1368 !
1369 !*  3.1       Preliminary calculations
1370 !
1371   ZCGAM    = COS(-ZRDSDG*XBETA)
1372   ZSGAM    = SIN(-ZRDSDG*XBETA)
1373   ZRACLAT0 = XRADIUS*COS(ZRDSDG*XLAT0)
1374   ZXE      = ZXBM0*ZCGAM+ZYBM0*ZSGAM            &
1375            - ZRACLAT0*(PLONOR-XLON0)*ZRDSDG  
1376   ZYE      =-ZXBM0*ZSGAM+ZYBM0*ZCGAM            &
1377            - ZRACLAT0*LOG(TAN(XPI/4.+PLATOR*ZRDSDG/2.))
1378 !
1379 !*  3.2       Conformal coordinates
1380 !
1381   ZXPR(:,:)   = ZRACLAT0*(ZLON(:,:)-XLON0)*ZRDSDG+ZXE
1382   ZYPR(:,:)   = ZRACLAT0*LOG(TAN(XPI/4.+PLAT(:,:)*ZRDSDG/2.))+ZYE
1383   !
1384   PXHATM(:,:) = ZXPR(:,:)*ZCGAM-ZYPR(:,:)*ZSGAM
1385   PYHATM(:,:) = ZXPR(:,:)*ZSGAM+ZYPR(:,:)*ZCGAM
1386 !
1387 !-------------------------------------------------------------------------------
1388 !
1389 !*  4.        EXIT
1390 !             ----
1391 !
1392 END IF
1393 RETURN
1394 !-------------------------------------------------------------------------------
1395 END SUBROUTINE SM_XYHAT_A
1396 !-------------------------------------------------------------------------------
1397 !
1398 !
1399 !-------------------------------------------------------------------------------
1400 !
1401 !*              6.   FUNCTION LATREF2
1402 !                    -----------------
1403 !-------------------------------------------------------------------------------
1404 !      #############################
1405        FUNCTION LATREF2(PLAT0,PRPK)
1406 !      #############################
1407 !
1408 !!****  *LATREF2 * - returns the Lambert second reference latitude
1409 !!
1410 !!     PURPOSE
1411 !!     -------
1412 !        This routine computes the second reference latitude 
1413 !      of a Lambert conformal projection for given projection
1414 !      parameter PRPK and primary reference latitude PLAT0.
1415 !        This second latitude is used in US and UK to define
1416 !      the secant Lambert projection (as a substitute for the
1417 !      cone constant PRPK used in France by IGN).
1418 !        This latitude is required to call the NCAR map projection
1419 !      package with the Lambert option.
1420 !
1421 !!**   METHOD
1422 !!     ------
1423 !!       The so-called "constant of the cone" equation is solved 
1424 !!     using a simple Newton-Raphson iteration. The spherical earth 
1425 !!     approximation is used.   
1426 !!
1427 !!       WARNING: ALL INPUT AND OUTPUT ANGLES ARE IN DEGREES...
1428 !!
1429 !!     EXTERNAL
1430 !!     --------
1431 !!       None
1432 !!
1433 !!     EXPLICIT ARGUMENTS
1434 !!     -------------------
1435 !!       PRPK    : projection factor       (no-unit, input)  
1436 !!       PLAT0   : map reference latitude  (degrees, input)             
1437 !!
1438 !!     IMPLICIT ARGUMENTS
1439 !!     ------------------
1440 !!       Module MODD_CST      : contains Physical constants
1441 !!          XPI        : Pi;    
1442 !!
1443 !!       Module MODD_LUNIT    : contains logical unit names
1444 !!          CLUOUT0    : Output listing file name
1445 !!
1446 !!     REFERENCE
1447 !!     ---------
1448 !!      Joly A., 1992, "Geographic parameters for ARPEGE/ALADIN",
1449 !!            Internal note ARPEGE/ALADIN, february 27,28p.
1450 !!      Levallois J., 1970, "Geodesie generale", Tome 2, Collection
1451 !!             de l'IGN, Eyrolles, Paris, 408p.
1452 !!      Pearson F. II, 1990,"Map projections: theory and applications",
1453 !!             CRC Press, Boca Raton, Florida, 372p. (Chapter 5).
1454 !!       
1455 !!     AUTHOR
1456 !!     ------
1457 !!      P.M.       *LA*
1458 !!
1459 !!     MODIFICATION
1460 !!     ------------
1461 !!       Original PM  24/05/94
1462 !!       Updated  PM  27/07/94
1463 !!       Updated  VD  25/08/94
1464 !!       Updated  VM  24/10/95 projection from north pole (XRPK<0)
1465 !!       Updated  VM  08/10/96 output-listing choice
1466 !!       Updated  IM  27/11/03 special case if projection plane is tangent
1467 !!
1468 !-------------------------------------------------------------------------------
1469 !
1470 !*     0.     DECLARATIONS
1471 !             ------------
1472 !
1473 USE MODD_CST
1474 USE MODD_LUNIT1
1475 !
1476 IMPLICIT NONE
1477 !
1478 !*     0.1    Declarations of arguments and results
1479 !
1480 REAL,INTENT(IN):: PLAT0,PRPK    ! Given first standard latitude (degrees)
1481                                 ! and projection parameter (cone 
1482                                 ! constant) for the Lambert conformal
1483                                 ! projection used.
1484 REAL :: LATREF2                 ! Returned latitude of the second 
1485                                 ! reference (or standard) parallel
1486                                 ! of the projection.
1487 !
1488 !*     0.2    Declarations of local variables
1489 !
1490 REAL    :: ZRPK
1491 REAL    :: ZRDSDG,ZEPSI,ZLAT0,ZLAT,ZDLAT,ZGLAT,ZGPRSG
1492 INTEGER :: ITER,ITERMAX
1493 INTEGER :: ILUOUT,IRESP  
1494 !
1495 !-------------------------------------------------------------------------------
1496 !
1497 !*     1.     PRELIMINARY CALCULATIONS
1498 !             ------------------------
1499 !
1500 ZRDSDG  = XPI/180.         ! Degree to radian conversion factor
1501 ZEPSI   = 10.*EPSILON(1.)   ! a small number
1502 ITERMAX = 10              ! number of iteration allowed
1503 !
1504 IF (PRPK ==SIN(ZLAT0*ZRDSDG)) THEN     ! projection plane tangent to the sphere
1505   LATREF2 = ZLAT0
1506 ELSE                                   !       "          intersect the sphere
1507 !
1508   ZLAT0   = PLAT0*ZRDSDG      ! Switch to radians
1509 !
1510   ZLAT    = XPI-4.*ATAN(SQRT((1.-PRPK)/(1.+PRPK)))-ZLAT0    
1511   ITER    = 0                  ! Choose the side of the nice root
1512   ZDLAT   = 0.                ! and sets up for the loop
1513 !
1514
1515 !
1516   IF (PRPK<0.) THEN    ! projection from north pole
1517     ZRPK=-PRPK
1518     ZLAT0=-ZLAT0
1519     ZLAT=-ZLAT
1520   ELSE                 ! projection from south pole
1521     ZRPK=PRPK
1522   ENDIF
1523 !
1524 !-------------------------------------------------------------------------------
1525 !
1526 !*     2.     NEWTON-RAPHSON LOOP
1527 !             -------------------
1528   DO
1529     ITER   = ITER+1
1530     ZLAT   = ZLAT+ZDLAT
1531     ZGLAT  =(COS(ZLAT)/COS(ZLAT0))*                         & 
1532             (((1.+SIN(ZLAT))/(1.+SIN(ZLAT0)))**(ZRPK/(1.-ZRPK)))
1533     ZGPRSG = ((ZRPK/(1.-ZRPK))*(COS(ZLAT)/(1.+SIN(ZLAT)))    &
1534             - (SIN(ZLAT)/COS(ZLAT)))*ZGLAT
1535     ZDLAT  = (1.-ZGLAT)/ZGPRSG
1536     !
1537     IF((ABS(ZGLAT-1.) <= ZEPSI).OR.(ITER >= ITERMAX))   EXIT
1538   END DO
1539 !
1540   IF (PRPK<0.) ZLAT=-ZLAT
1541   LATREF2  = ZLAT/ZRDSDG     ! Degrees restored
1542 !
1543 ENDIF
1544 !-------------------------------------------------------------------------------
1545 !
1546 !*  3.        EXIT
1547 !             ----
1548 !
1549 IF(ITER <= ITERMAX)  RETURN
1550 !
1551 CALL FMLOOK(CLUOUT,CLUOUT,ILUOUT,IRESP)
1552 !
1553 WRITE(ILUOUT,*) ' Error in function LATREF2 (module MODE_GRIDPROJ)'
1554 WRITE(ILUOUT,*) ' Function fails to converge after ',ITER,' iterations.'
1555 WRITE(ILUOUT,*) ' LATREF2=',LATREF2,' Residual=',ZGLAT-1.,           &
1556                 ' ZEPSI=',ZEPSI,' Last increment=',ZDLAT/ZRDSDG
1557 WRITE(ILUOUT,*) ' JOB ABORTS...'
1558 STOP
1559 !-------------------------------------------------------------------------------
1560 END FUNCTION LATREF2
1561 !-------------------------------------------------------------------------------
1562 !
1563 END MODULE MODE_GRIDPROJ