SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
mode_gridtype_gauss.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_gauss
7 !##########################
8 !
9 !############################################################################
10 !############################################################################
11 !############################################################################
12 !
13 USE yomhook ,ONLY : lhook, dr_hook
14 USE parkind1 ,ONLY : jprb
15 !
16 USE modi_abor1_sfx
17 !
18 IMPLICIT NONE
19 !
20  CONTAINS
21 !
22 !############################################################################
23 !############################################################################
24 !############################################################################
25 ! ####################################################################
26  SUBROUTINE put_gridtype_gauss(PGRID_PAR,KNLATI, PLAPO,PLOPO,PCODIL,KNLOPA, &
27  kl,plat,plon,plat_xy,plon_xy,pmesh_size , &
28  ploninf,platinf,plonsup,platsup )
29 ! ####################################################################
30 !
31 !!**** *PUT_GRIDTYPE_GAUSS* - routine to store in PGRID_PAR the horizontal grid
32 !!
33 !! AUTHOR
34 !! ------
35 !! V. Masson *Meteo France*
36 !!
37 !! MODIFICATIONS
38 !! -------------
39 !! Original 01/2004
40 !! (B. Decharme) 2008 multiples changes
41 !-------------------------------------------------------------------------------
42 !
43 !* 0. DECLARATIONS
44 ! ------------
45 !
46 IMPLICIT NONE
47 !
48 !
49 !* 0.1 Declarations of arguments
50 ! -------------------------
51 !
52 INTEGER, INTENT(IN) :: knlati ! number of pseudo-latitudes
53 REAL, INTENT(IN) :: plapo ! latitude of the rotated pole (deg)
54 REAL, INTENT(IN) :: plopo ! logitude of the rotated pole (rad)
55 REAL, INTENT(IN) :: pcodil ! stretching factor
56 INTEGER, DIMENSION(KNLATI), INTENT(IN) :: knlopa ! number of pseudo-longitudes
57 ! ! on each pseudo-latitude circle
58 ! ! on pseudo-northern hemisphere
59 ! ! (starting from the rotated pole)
60 INTEGER, INTENT(IN) :: kl ! number of points used
61 REAL, DIMENSION(:), INTENT(IN) :: plat ! latitudes of points
62 REAL, DIMENSION(:), INTENT(IN) :: plon ! longitudes of points
63 REAL, DIMENSION(:), INTENT(IN) :: plat_xy ! pseudo-latitudes of points
64 REAL, DIMENSION(:), INTENT(IN) :: plon_xy ! pseudo-longitudes of points
65 REAL, DIMENSION(:), INTENT(IN) :: pmesh_size ! Mesh size
66 ! _____ Sup
67 REAL, DIMENSION(:), INTENT(IN) :: platsup ! Grid corner Latitude | |
68 REAL, DIMENSION(:), INTENT(IN) :: plonsup ! Grid corner Longitude | |
69 REAL, DIMENSION(:), INTENT(IN) :: platinf ! Grid corner Latitude |_____|
70 REAL, DIMENSION(:), INTENT(IN) :: ploninf ! Grid corner Longitude Inf
71 !
72 REAL, DIMENSION(:), POINTER :: pgrid_par ! parameters defining this grid
73 REAL(KIND=JPRB) :: zhook_handle
74 !
75 !
76 !
77 !* 0.2 Declarations of local variables
78 ! -------------------------------
79 !
80 !-------------------------------------------------------------------------------
81 !
82 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:PUT_GRIDTYPE_GAUSS',0,zhook_handle)
83 ALLOCATE(pgrid_par(5+knlati+9*kl))
84 pgrid_par(1) = knlati
85 pgrid_par(2) = plapo
86 pgrid_par(3) = plopo
87 pgrid_par(4) = pcodil
88 pgrid_par(5:4+knlati)= knlopa(:)
89 pgrid_par(5+knlati) = kl
90 pgrid_par(6+knlati:5+knlati+kl) = plat(:)
91 pgrid_par(6+knlati+kl:5+knlati+2*kl) = plon(:)
92 pgrid_par(6+knlati+2*kl:5+knlati+3*kl) = plat_xy(:)
93 pgrid_par(6+knlati+3*kl:5+knlati+4*kl) = plon_xy(:)
94 pgrid_par(6+knlati+4*kl:5+knlati+5*kl) = pmesh_size(:)
95 pgrid_par(6+knlati+5*kl:5+knlati+6*kl) = ploninf(:)
96 pgrid_par(6+knlati+6*kl:5+knlati+7*kl) = platinf(:)
97 pgrid_par(6+knlati+7*kl:5+knlati+8*kl) = plonsup(:)
98 pgrid_par(6+knlati+8*kl:5+knlati+9*kl) = platsup(:)
99 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:PUT_GRIDTYPE_GAUSS',1,zhook_handle)
100 !
101 !-------------------------------------------------------------------------------
102 END SUBROUTINE put_gridtype_gauss
103 !############################################################################
104 !############################################################################
105 !############################################################################
106 ! ####################################################################
107  SUBROUTINE get_gridtype_gauss(PGRID_PAR,KNLATI, &
108  plapo,plopo,pcodil,knlopa,kl, &
109  plat,plon,plat_xy,plon_xy,pmesh_size,&
110  ploninf,platinf,plonsup,platsup )
111 ! ####################################################################
112 !
113 !!**** *GET_GRIDTYPE_GAUSS* - routine to get from PGRID_PAR the horizontal grid
114 !!
115 !! AUTHOR
116 !! ------
117 !! V. Masson *Meteo France*
118 !!
119 !! MODIFICATIONS
120 !! -------------
121 !! Original 01/2004
122 !-------------------------------------------------------------------------------
123 !
124 !* 0. DECLARATIONS
125 ! ------------
126 !
127 IMPLICIT NONE
128 !
129 !
130 !* 0.1 Declarations of arguments
131 ! -------------------------
132 !
133 REAL, DIMENSION(:), INTENT(IN) :: pgrid_par! parameters defining this grid
134 INTEGER, INTENT(OUT), OPTIONAL :: knlati ! number of pseudo-latitudes
135 REAL, INTENT(OUT), OPTIONAL :: plapo ! latitude of the rotated pole (deg)
136 REAL, INTENT(OUT), OPTIONAL :: plopo ! logitude of the rotated pole (deg)
137 REAL, INTENT(OUT), OPTIONAL :: pcodil ! stretching factor
138 INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: knlopa ! number of pseudo-longitudes
139 ! ! on each pseudo-latitude circle
140 ! ! on pseudo-northern hemisphere
141 ! ! (starting from the rotated pole)
142 INTEGER, INTENT(OUT), OPTIONAL :: kl ! number of points
143 REAL, DIMENSION(:), OPTIONAL, INTENT(OUT) :: plat ! latitude
144 REAL, DIMENSION(:), OPTIONAL, INTENT(OUT) :: plon ! longitude
145 REAL, DIMENSION(:), OPTIONAL, INTENT(OUT) :: plat_xy ! pseudo-latitude
146 REAL, DIMENSION(:), OPTIONAL, INTENT(OUT) :: plon_xy ! pseudo-longitude
147 REAL, DIMENSION(:), OPTIONAL, INTENT(OUT) :: pmesh_size ! Mesh size
148 ! _____ Sup
149 REAL, DIMENSION(:), OPTIONAL, INTENT(OUT) :: platsup ! Grid corner Latitude | |
150 REAL, DIMENSION(:), OPTIONAL, INTENT(OUT) :: plonsup ! Grid corner Longitude | |
151 REAL, DIMENSION(:), OPTIONAL, INTENT(OUT) :: platinf ! Grid corner Latitude |_____|
152 REAL, DIMENSION(:), OPTIONAL, INTENT(OUT) :: ploninf ! Grid corner Longitude Inf
153 !
154 !
155 !* 0.2 Declarations of local variables
156 ! -------------------------------
157 !
158 INTEGER :: il
159 INTEGER :: inlati
160 REAL(KIND=JPRB) :: zhook_handle
161 !-------------------------------------------------------------------------------
162 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:GET_GRIDTYPE_GAUSS',0,zhook_handle)
163 IF (present(knlati)) knlati = pgrid_par(1)
164 IF (present(plapo)) plapo = pgrid_par(2)
165 IF (present(plopo)) plopo = pgrid_par(3)
166 IF (present(pcodil)) pcodil = pgrid_par(4)
167 !
168 inlati = pgrid_par(1)
169 !
170 IF (present(knlopa)) THEN
171  knlopa(:) = pgrid_par(5:4+inlati)
172 END IF
173 !
174 il = pgrid_par(5+inlati)
175 !
176 IF (present(kl)) THEN
177  kl = il
178 END IF
179 !
180 IF (present(plat)) THEN
181  IF (SIZE(plat)/=il) THEN
182  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: WRONG SIZE FOR PLAT')
183  END IF
184  plat(:) = pgrid_par(6+inlati:5+inlati+il)
185 END IF
186 !
187 IF (present(plon)) THEN
188  IF (SIZE(plon)/=il) THEN
189  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: WRONG SIZE FOR PLON')
190  END IF
191  plon(:) = pgrid_par(6+inlati+il:5+inlati+2*il)
192 END IF
193 !
194 IF (present(plat_xy)) THEN
195  IF (SIZE(plat_xy)/=il) THEN
196  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: WRONG SIZE FOR PLAT_XY')
197  END IF
198  plat_xy(:) = pgrid_par(6+inlati+2*il:5+inlati+3*il)
199 END IF
200 !
201 IF (present(plon_xy)) THEN
202  IF (SIZE(plon_xy)/=il) THEN
203  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: WRONG SIZE FOR PLON_XY')
204  END IF
205  plon_xy(:) = pgrid_par(6+inlati+3*il:5+inlati+4*il)
206 END IF
207 !
208 IF (present(pmesh_size)) THEN
209  IF (SIZE(pmesh_size)/=il) THEN
210  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: WRONG SIZE FOR PMESH_SIZE')
211  END IF
212  pmesh_size(:) = pgrid_par(6+inlati+4*il:5+inlati+5*il)
213 END IF
214 !
215 IF (present(ploninf)) THEN
216  IF (SIZE(ploninf)/=il) THEN
217  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: WRONG SIZE FOR PLONINF')
218  END IF
219  ploninf(:) = pgrid_par(6+inlati+5*il:5+inlati+6*il)
220 END IF
221 !
222 IF (present(platinf)) THEN
223  IF (SIZE(platinf)/=il) THEN
224  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: WRONG SIZE FOR PLATINF')
225  END IF
226  platinf(:) = pgrid_par(6+inlati+6*il:5+inlati+7*il)
227 END IF
228 !
229 IF (present(plonsup)) THEN
230  IF (SIZE(plonsup)/=il) THEN
231  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: WRONG SIZE FOR PLONSUP')
232  END IF
233  plonsup(:) = pgrid_par(6+inlati+7*il:5+inlati+8*il)
234 END IF
235 !
236 IF (present(platsup)) THEN
237  IF (SIZE(platsup)/=il) THEN
238  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: WRONG SIZE FOR PLATSUP')
239  END IF
240  platsup(:) = pgrid_par(6+inlati+8*il:5+inlati+9*il)
241 END IF
242 !
243 !
244 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:GET_GRIDTYPE_GAUSS',1,zhook_handle)
245 !
246 !-------------------------------------------------------------------------------
247 END SUBROUTINE get_gridtype_gauss
248 !############################################################################
249 !############################################################################
250 !############################################################################
251 ! ####################################################################
252  SUBROUTINE latlon_gauss(PLON_XY,PLAT_XY,KL,PLOPO,PLAPO,PCODIL,PLON,PLAT)
253 ! ####################################################################
254 !
255 !!**** *LATLON_GAUSS* - computes the coordinates on the real sphere
256 !!
257 !! AUTHOR
258 !! ------
259 !! F. Taillefer *Meteo France*
260 !!
261 !! MODIFICATIONS
262 !! -------------
263 !! Original 10/2007
264 !-------------------------------------------------------------------------------
265 !
266 !* 0. DECLARATIONS
267 ! ------------
268 !
269 USE modd_csts, ONLY : xpi
270 !
271 IMPLICIT NONE
272 !
273 !
274 !* 0.1 Declarations of arguments
275 ! -------------------------
276 !
277 INTEGER, INTENT(IN) :: kl ! number of grid points
278 REAL, INTENT(IN) :: plapo ! pole latitude
279 REAL, INTENT(IN) :: plopo ! pole longitude
280 REAL, INTENT(IN) :: pcodil ! stretching factor
281 REAL, DIMENSION(KL), INTENT(IN) :: plat_xy ! pseudo-latitudes of points (deg)
282 REAL, DIMENSION(KL), INTENT(IN) :: plon_xy ! pseudo-longitudes of points (deg)
283 REAL, DIMENSION(KL), INTENT(OUT):: plat ! latitudes of points (deg)
284 REAL, DIMENSION(KL), INTENT(OUT):: plon ! longitudes of points (deg)
285 !
286 !
287 !* 0.2 Declarations of local variables
288 ! -------------------------------
289 ! called 1 variables : on the transformed sphere (rotated ans stretched)
290 ! called 2 variables : on the rotated but no stretched sphere
291 ! called 3 variables : on the real sphere
292 
293 INTEGER :: jp
294 REAL :: zclo3,zconr,zinterm,zlat1,zlat2,zlat3,zlon1,zlon2,zlon3
295 REAL :: zlatp,zlonp,zsla3,zslo3,zr
296 REAL(KIND=JPRB) :: zhook_handle
297 !-------------------------------------------------------------------------------
298 !
299 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:LATLON_GAUSS',0,zhook_handle)
300 !
301 zconr=xpi/180.
302 !
303 zlonp=zconr*plopo
304 zlatp=zconr*plapo
305 !
306 DO jp = 1,kl
307  zlon1=zconr*plon_xy(jp)
308  zlat1=zconr*plat_xy(jp)
309 ! move from the stretched to the no stretched sphere
310  zinterm=1./pcodil*cos(zlat1)/max(1.e-12,1.+sin(zlat1))
311  zlat2=2.*atan((1.-zinterm)/(1.+zinterm))
312  zlon2=zlon1
313 ! move from the rotated sphere to the real one
314 !
315 ! calculation of the latitude
316  zsla3=-cos(zlat2)*cos(zlon2)*cos(zlatp)+sin(zlat2)*sin(zlatp)
317  IF (zsla3>1. .OR. zsla3<-1.) THEN
318  WRITE(0,*) 'be carefull --> sinus >1'
319  zsla3=min(1.,max(-1.,zsla3))
320  ENDIF
321  zlat3=asin(zsla3)
322 !
323 ! calculation of the sine and cosine of the longitude
324  zslo3=(cos(zlat2)*cos(zlon2)*sin(zlatp)*sin(zlonp)&
325  +cos(zlat2)*sin(zlon2)*cos(zlonp)&
326  +sin(zlat2)*cos(zlatp)*sin(zlonp)) / cos(zlat3)
327  zclo3=(cos(zlat2)*cos(zlon2)*sin(zlatp)*cos(zlonp)&
328  -cos(zlat2)*sin(zlon2)*sin(zlonp)&
329  +sin(zlat2)*cos(zlatp)*cos(zlonp)) / cos(zlat3)
330 !
331 ! Conversion from rectangular to polar to get the longitude
332  zr=sqrt(zclo3*zclo3+zslo3*zslo3)
333  IF (zclo3==0.) THEN
334  IF (zslo3==0.) THEN
335  zlon3=0.
336  ELSEIF (zslo3>0.) THEN
337  zlon3=xpi/2.
338  ELSE
339  zlon3=-xpi/2.
340  ENDIF
341  ELSE
342  zinterm=atan(zslo3/zclo3)
343  IF (zclo3>=0.) THEN
344  zlon3=zinterm
345  ELSEIF (zslo3>=0.) THEN
346  zlon3=zinterm+xpi
347  ELSE
348  zlon3=zinterm-xpi
349  ENDIF
350  ENDIF
351 !
352 ! Conversion from radians to degrees
353  plon(jp)=zlon3/zconr
354  IF (plon(jp)<0.) plon(jp)=plon(jp)+360.
355  plat(jp)=zlat3/zconr
356 !
357 END DO
358 !
359 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:LATLON_GAUSS',1,zhook_handle)
360 !
361 !-------------------------------------------------------------------------------
362 END SUBROUTINE latlon_gauss
363 !############################################################################
364 !############################################################################
365 !############################################################################
366 ! ####################################################################
367  SUBROUTINE comp_gridtype_gauss(KNLATI,KNLOPA,KL,KTYP,PLAT_XY,PLON_XY)
368 ! ####################################################################
369 !
370 !!**** *COMP_GRIDTYPE_GAUSS* - computes the gaussian grid
371 !!
372 !! AUTHOR
373 !! ------
374 !! V. Masson *Meteo France*
375 !!
376 !! MODIFICATIONS
377 !! -------------
378 !! Original 10/2006
379 !! F.Taillefer 10/2007 : adapt gaussian latitudes calculation
380 !-------------------------------------------------------------------------------
381 !
382 !* 0. DECLARATIONS
383 ! ------------
384 !
385 USE eggangles , ONLY : p_asin
386 USE modd_csts, ONLY : xpi
387 !
388 IMPLICIT NONE
389 !
390 !
391 !* 0.1 Declarations of arguments
392 ! -------------------------
393 !
394 INTEGER, INTENT(IN) :: knlati ! number of pseudo-latitudes
395 INTEGER, DIMENSION(KNLATI), INTENT(IN) :: knlopa ! number of pseudo-longitudes
396 ! ! on each pseudo-latitude circle
397 ! ! on pseudo-northern hemisphere
398 ! ! (starting from the rotated pole)
399 INTEGER, INTENT(IN) :: kl ! number of points used
400 INTEGER, INTENT(IN) :: ktyp ! type of transform
401 REAL, DIMENSION(KL), INTENT(INOUT) :: plat_xy ! pseudo-latitudes of points (deg)
402 REAL, DIMENSION(KL), INTENT(INOUT) :: plon_xy ! pseudo-longitudes of points (deg)
403 !
404 !
405 !* 0.2 Declarations of local variables
406 ! -------------------------------
407 !
408 INTEGER :: jl ! point loop counter
409 INTEGER :: jx ! longitude loop counter
410 INTEGER :: jy ! latitude loop counter
411 REAL :: zrd, zi
412 REAL, DIMENSION(KNLATI) :: znlopa, zsinla, zwg
413 REAL, DIMENSION(KNLATI) :: zdsinla
414 REAL(KIND=JPRB) :: zhook_handle
415 !-------------------------------------------------------------------------------
416 !
417 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:COMP_GRIDTYPE_GAUSS',0,zhook_handle)
418 !
419 zrd = 180. / xpi
420 zi=0.
421 IF (ktyp==1) zi=180.
422 !
423 znlopa=float(knlopa)
424 !
425 ! gaussian latitudes calculation
426  CALL latitudes_gauss(knlati,zsinla,zdsinla,zwg)
427 !
428 jl=knlati/2
429 DO jy=1,jl
430  zwg(jy)= p_asin(zsinla(jy))
431  zwg(knlati+1-jy)= -p_asin(zsinla(jy))
432 END DO
433 !
434 jl=0
435 !* loop on latitudes (from north pole)
436 DO jy=1,knlati
437 !* loop on longitudes (from 0°, to the east)
438  DO jx=1,knlopa(jy)
439  jl=jl+1
440  plat_xy(jl) = zwg(jy)*zrd
441  plon_xy(jl) = 360. * float(jx-1) / znlopa(jy) + zi
442  END DO
443 END DO
444 !
445 IF (jl/=kl) THEN
446  WRITE(0,*) ' PB in the total number of points of the gaussian grid '
447  WRITE(0,*) ' check your namelist and rerun !'
448  CALL abor1_sfx('MODE_GRIDTYPE_GAUSS: PB IN THE TOTAL NUMBER OF POINTS')
449 ENDIF
450 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:COMP_GRIDTYPE_GAUSS',1,zhook_handle)
451 !
452 !-------------------------------------------------------------------------------
453 END SUBROUTINE comp_gridtype_gauss
454 !############################################################################
455 !############################################################################
456 !############################################################################
457 ! ####################################################################
458  SUBROUTINE gauss_grid_limits(KNLATI,KNLOPA,PXINF,PXSUP,PYINF,PYSUP)
459 ! ####################################################################
460 !
461 !!**** *GAUSS_GRID_LIMITS* - computes the gaussian grid "boxes"
462 !!
463 !! AUTHOR
464 !! ------
465 !! V. Masson *Meteo France*
466 !!
467 !! MODIFICATIONS
468 !! -------------
469 !! Original 10/2006
470 !! F. Taillefer 08/2007 pb with the gaussian latitudes
471 !-------------------------------------------------------------------------------
472 !
473 !* 0. DECLARATIONS
474 ! ------------
475 !
476 USE eggangles , ONLY : p_asin
477 USE modd_csts, ONLY : xpi
478 !
479 IMPLICIT NONE
480 !
481 !* 0.1 Declarations of arguments
482 ! -------------------------
483 !
484 INTEGER, INTENT(IN) :: knlati! number of pseudo-latitudes
485 INTEGER, DIMENSION(KNLATI),INTENT(IN) :: knlopa! number of pseudo-longitudes
486 ! ! on each pseudo-latitude circle
487 ! ! on pseudo-northern hemisphere
488 ! ! (starting from the rotated pole)
489 REAL, DIMENSION(:), INTENT(OUT):: pxinf ! minimum pseudo longitude of the grid point (deg)
490 REAL, DIMENSION(:), INTENT(OUT):: pxsup ! maximum pseudo longitude of the grid point (deg)
491 REAL, DIMENSION(:), INTENT(OUT):: pyinf ! minimum pseudo latitude of the grid point (deg)
492 REAL, DIMENSION(:), INTENT(OUT):: pysup ! maximum pseudo latitude of the grid point (deg)
493 !
494 !
495 !* 0.2 Declarations of local variables
496 ! -------------------------------
497 !
498 INTEGER :: jl ! point loop counter
499 INTEGER :: jx ! longitude loop counter
500 INTEGER :: jy ! latitude loop counter
501 REAL :: znlati, zrd
502 REAL, DIMENSION(KNLATI) :: znlopa, zsinla, zwg
503 REAL, DIMENSION(KNLATI) :: zdsinla
504 REAL(KIND=JPRB) :: zhook_handle
505 !-------------------------------------------------------------------------------
506 !
507 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:GAUSS_GRID_LIMITS',0,zhook_handle)
508 !
509 zrd = 180. / xpi
510 !
511 znlati=float(knlati)
512 znlopa=float(knlopa)
513 !
514 ! gaussian latitudes calculation
515  CALL latitudes_gauss(knlati,zsinla,zdsinla,zwg)
516 !
517 jl=knlati/2
518 DO jy=1,jl
519  zwg(jy)= p_asin(zsinla(jy))*zrd
520  zwg(knlati+1-jy)= -p_asin(zsinla(jy))*zrd
521 END DO
522 !
523 jl=0
524 !
525 !* loop on latitudes (from north pole)
526 DO jy=1,knlati
527 !* loop on longitudes (from 0°, to the east)
528  DO jx=1,knlopa(jy)
529  jl=jl+1
530  IF (jy==1) THEN
531  pysup(jl) = 90.
532  ELSE
533  pysup(jl)= zwg(jy)+(zwg(jy-1)-zwg(jy))/2.
534  ENDIF
535  IF (jy==knlati) THEN
536  pyinf(jl) = -90.
537  ELSE
538  pyinf(jl) = zwg(jy)-(zwg(jy)-zwg(jy+1))/2.
539  ENDIF
540  pxsup(jl) = 360. * (float(jx)-0.5) / znlopa(jy)
541  pxinf(jl) = 360. * (float(jx)-1.5) / znlopa(jy)
542  END DO
543 END DO
544 IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:GAUSS_GRID_LIMITS',1,zhook_handle)
545 !
546 !-------------------------------------------------------------------------------
547 END SUBROUTINE gauss_grid_limits
548 !############################################################################
549  !############################################################################
550  ! ####################################################################
551 SUBROUTINE xy_gauss(PCODIL,KSIZE,PNODATA,PVALUE,PLAT_XY,PLON_XY)
552  ! ####################################################################
553  !
554  !!**** *LATLON_GAUSS * - Routine to compute coordinates on a transform sphere
555  !! from geographical coordinates
556  !!
557  !! PURPOSE
558  !! -------
559  ! This routine computes the latitude and longitude a real coordinates
560  ! given array to an arpege model coordinates (rotated stretched)
561  !
562  !
563  !
564  !!** METHOD
565  !! ------
566  !! use of rotations routines (eggmrt) and streching conformal formulae
567  !! to pass from real sphere (PLAT,PLON) transform sphere (PLAT_XY, PLON_XY)
568  !!
569  !! EXTERNAL
570  !! --------
571  !! None
572  !!
573  !! REFERENCE
574  !! ---------
575  !! Arpege DOC "Sphere Transphormee" Chapitre 7 version du 4/6/1991
576  !! J-D Gril for GEO_GAUSS 2005
577  !! J-D Gril Doc for EGGANGLES routines (new EGGX) 2005
578  !!
579  !! AUTHOR
580  !! ------
581  !! J-D Gril
582  !!
583  !! MODIFICATION
584  !! ------------
585  !! Original 10/2005
586  !
587  !-------------------------------------------------------------------------------
588  !
589  !* 0. DECLARATIONS
590  ! ------------
591  !
592  USE modd_get_mesh_index_gauss, ONLY : xlon, xlat, xcost, xsintc, xsints, xcosn, xsinn, &
593  xlonp, xlatp, xcosp, xsinp, xpi, x1, x2, xdr
594  !
595  IMPLICIT NONE
596  !
597  !* 0.1 Declarations of arguments and results
598  !
599  REAL, INTENT(IN) :: pcodil
600  INTEGER, INTENT(IN) :: ksize
601  REAL, INTENT(IN) :: pnodata
602  REAL, DIMENSION(:), INTENT(IN) :: pvalue ! value of the point to add
603  REAL, DIMENSION(:), INTENT(OUT):: plat_xy,plon_xy
604  !
605  !* 0.2 Declarations of local variables
606  !
607  REAL :: zcos1, zv1, zinvs, zlat1, zcos2, zlon1, zinvc, zsinn
608  REAL :: zv2, zsint, zcost, zv3, zlat2, zm, zlon2
609  !
610  INTEGER :: jj, idn, idt
611  !
612  REAL(KIND=JPRB) :: zhook_handle
613  !--------------------------------------------------------------------------------
614  !
615  !* 1. Preliminary calculations
616  ! ------------------------
617  !
618  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:XY_GAUSS',0,zhook_handle)
619  !
620 !
621 !$OMP PARALLEL DO PRIVATE(IDN,IDT,ZCOS1, ZV1, ZLAT1, ZCOS2, ZLON1, ZINVC, &
622 !$OMP ZSINN, ZV2, ZINVS, ZSINT, ZCOST, ZV3, ZLAT2, ZM, ZLON2)
623  DO jj=1,SIZE(pvalue)
624  !
625  IF (pvalue(jj)==pnodata) cycle
626  !
627  idn = mod(jj,ksize)
628  IF (idn==0) idn=ksize
629  idt = ceiling(1.*jj/ksize)
630  !
631  zcos1 = xcosn(idn)*xcost(idt)
632  !ZCOS1 = COS(XLAT(IDT))*COS(XLON(IDN)-XLONP)
633  zv1 = min(1.,max(-1.,xsints(idt)+xcosp*zcos1))
634  !ZV1 = MIN(1.,MAX(-1.,SIN(XLATP)*SIN(XLAT(IDT))+COS(XLATP)*ZCOS1))
635  zlat1 = asin(zv1)
636  !
637  zcos2 = cos(zlat1)
638  !
639  zlon1 = 0.0
640  IF (zcos2 /= 0.0) THEN
641  !ZINVC = 1./ZCOS2
642  !ZSINN = - XCOST(IDT)*XSINN(IDN)*ZINVC
643  !ZV2 = MIN(1.,MAX(-1.,(XSINTC(IDT)-XSINP*ZCOS1)*ZINVC))
644  zsinn = - xcost(idt)*xsinn(idn)/zcos2
645  !ZSINN = - COS(XLAT(IDT))*SIN(XLON(IDN)-XLONP)/ZCOS2
646  zv2 = min(1.,max(-1.,(xsintc(idt)-xsinp*zcos1)/zcos2))
647  !ZV2 = MIN(1.,MAX(-1.,(COS(XLATP)*SIN(XLAT(IDT))-SIN(XLATP)*ZCOS1)/ZCOS2))
648  zlon1 = acos(zv2) * sign(1.,zsinn)
649  ENDIF
650  !
651  !-----------------------
652  !ZINVS = 1./(X2+X1*ZV1)
653  !ZSINT = ZINVS*(X1+X2*ZV1)
654  !ZCOST = ZINVS*ZCOS2*PCODIL*2.0
655  zsint = (x1+x2*zv1)/(x2+x1*zv1)
656  zcost = 2.0*pcodil*zcos2/(x2+x1*zv1)
657  !
658  zv3 = min(1.,max(-1.,zcost))
659  zlat2 = acos(zv3)*sign(1.,zsint)
660  !
661  zm = mod(zlon1,xpi)
662  zlon2 = (zm-xpi*mod(REAL(INT(ZLON1/XPI)),2.0))*sign(1.0,zlon1)*sign(1.0,zm)
663  !
664  !---------------------------------------------------------------------------------
665  !
666  !* 3. EXIT
667  ! ----
668  !
669  plat_xy(jj) = zlat2 / xdr
670  plon_xy(jj) = zlon2 / xdr
671  !
672  IF (abs(plon_xy(jj)-360.)<1.e-4) plon_xy(jj) = 0.
673  !
674  ENDDO
675 !$OMP END PARALLEL DO
676 !
677  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:XY_GAUSS',1,zhook_handle)
678  !---------------------------------------------------------------------------------
679  END SUBROUTINE xy_gauss
680  !---------------------------------------------------------------------------------
681  !############################################################################
682  ! ####################################################################
683  SUBROUTINE map_factor_gauss(PLAPO,PLOPO,PCODIL, &
684  plat,plon,pmap)
685  ! ####################################################################
686  !
687  !!**** *MAP_FACTOR_GAUSS * - Routine to compute map factor for points
688  !! in real sphere coordinates
689  !!
690  !!
691  !! PURPOSE
692  !! -------
693  !
694  !! REFERENCE
695  !! ---------
696  !! Arpege DOC "Sphere Transphormee" Chapitre 7 version du 4/6/1991
697  !! J-D Gril for GEO_GAUSS 2005
698  !! J-D Gril Doc for EGGANGLES routines (new EGGX) 2005
699  !!
700  !! AUTHOR
701  !! ------
702  !! J-D Gril
703  !!
704  !! MODIFICATION
705  !! ------------
706  !! Original 10/2005
707  !!
708  !-------------------------------------------------------------------------------
709  !
710  !* 0. DECLARATIONS
711  ! ------------
712  !
713  USE eggangles,ONLY : lola, angle_domain
714  USE mode_geo_gauss,ONLY : map_fac
715  !
716  IMPLICIT NONE
717  !
718  !* 0.1 Declarations of arguments and results
719  !
720  REAL, INTENT(IN) :: plapo ! latitude of pole
721  REAL, INTENT(IN) :: plopo ! longitude of pole
722  REAL, INTENT(IN) :: pcodil ! coefficient of dilatation
723  REAL, DIMENSION(:), INTENT(IN) :: plat ! given geographic latitudes
724  REAL, DIMENSION(:), INTENT(IN) :: plon ! given geographic longitudes
725  !
726  REAL, DIMENSION(SIZE(PLAT)), INTENT(OUT):: pmap ! map factor
727  !
728  !* 0.2 Declarations of local variables
729  !
730  !
731  TYPE(lola) :: tzpole
732  TYPE(lola), DIMENSION(SIZE(PLAT)) :: tzptci
733  !
734  REAL :: zpi
735  REAL :: zdr
736  REAL(KIND=JPRB) :: zhook_handle
737  !-------------------------------------------------------------------------------
738  !
739  !* 1. PRELIMINARY CALCULATION
740  ! -----------------------
741  !
742  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:MAP_FACTOR_GAUSS',0,zhook_handle)
743  zpi = 4.*atan(1.)
744  zdr = zpi / 180.
745  !
746  tzptci(:)%LON = angle_domain(plon(:),dom='0+',unit='D') * zdr
747  tzptci(:)%LAT = plat(:) * zdr
748  tzpole%LON = angle_domain(plopo,dom='0+',unit='D') * zdr
749  tzpole%LAT = plapo * zdr
750  !
751  !-------------------------------------------------------------------------------
752  !
753  !* 2. Calcul
754  ! ------
755  !
756  pmap(:) = map_fac(tzpole,pcodil,tzptci)
757  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:MAP_FACTOR_GAUSS',1,zhook_handle)
758  !
759  !-------------------------------------------------------------------------------
760  END SUBROUTINE map_factor_gauss
761  !-------------------------------------------------------------------------------
762  !############################################################################
763  !##################################
764  SUBROUTINE latitudes_gauss(KN,PL,PDL,PW)
765  !##################################
766 
767  !**** *SUGAW36 * - Routine to initialize the Gaussian
768  ! abcissa and the associated weights
769  ! Frozen CY36 version.
770 
771  ! Purpose.
772  ! --------
773  ! Initialize arrays PL,PDL and PW (quadrature abscissas and weights)
774  !** Interface.
775  ! ----------
776  ! *CALL* *SUGAW36(...) *
777 
778  ! Explicit arguments :
779  ! --------------------
780  ! INPUT:
781  ! KN : Number of Gauss abscissas
782 
783  ! OUTPUT:
784  ! PL (KN) : Abscissas of Gauss
785  ! PDL(KN) : Idem in quadruple precision (OPTIONNAL)
786  ! PW (KN) : Weights of the Gaussian integration
787 
788  ! PL (i) is the abscissa i starting from the northern pole, it is
789  ! the cosine of the colatitude of the corresponding row of the collocation grid.
790 
791  ! KOUTENV : Contains printings environment variables.
792  ! KOUTENV(1): output logical unit.
793  ! KOUTENV(2): option for level of printings.
794  !
795  ! Implicit arguments :
796  ! --------------------
797  ! None
798 
799  ! Method.
800  ! -------
801  ! See documentation
802 
803  ! Externals.
804  ! ----------
805 
806  ! Reference.
807  ! ----------
808 
809  ! S.L. Belousov, Tables of normalized associated Legendre Polynomials, Pergamon Press (1962)
810  ! P.N. Swarztrauber, On computing the points and weights for Gauss-Legendre quadrature,
811  ! SIAM J. Sci. Comput. Vol. 24 (3) pp. 945-954 (2002)
812 
813  ! Author.
814  ! -------
815  ! Mats Hamrud and Philippe Courtier *ECMWF*
816  ! Original : 87-10-15
817 
818  ! Modifications.
819  ! --------------
820  ! Philippe Courtier : 92-12-19 Multitasking
821  ! Ryad El Khatib : 94-04-20 Remove unused comdecks pardim and yomdim
822  ! Mats Hamrud : 94-08-12 Printing level
823  ! K. Yessad (Sep 2008): cleaning, improve comments.
824  ! K. YESSAD (NOV 2008): make consistent arp/SUGAW36 and tfl/SUGAW.
825  ! Nils Wedi + Mats Hamrud, 2009-02-05 revised following Swarztrauber, 2002
826  ! K. Yessad (June 2009): externalisation + in-lining.
827  ! GCO + K. Yessad (Dec 2012): restore quadruple precision features.
828  ! K. Yessad (Dec 2012): simplify use of SUGAW36.
829  ! ------------------------------------------------------------------
830 
831 
832  ! ------------------------------------------------------------------
833 
834  IMPLICIT NONE
835 
836  INTEGER,INTENT(IN) :: kn
837  REAL,INTENT(OUT) :: pl(kn)
838  REAL,INTENT(OUT) :: pdl(kn)
839  REAL,INTENT(OUT) :: pw(kn)
840 
841  ! ------------------------------------------------------------------
842 
843  REAL, DIMENSION(0:KN,0:KN) :: zfn
844  REAL, DIMENSION(0:KN/2) :: zfnlat
845  REAL, DIMENSION(KN) :: zreg, zli, zt, zmod, zm, zr
846 
847  REAL :: zfnn, zacos
848  REAL :: zpi, z, zra, zeps, zw, zx, zxn
849  REAL :: zdlx, zdlxn, zdlk, zdlldn, zzdlxn, zdlmod
850 
851  INTEGER, DIMENSION(KN) :: iter
852 
853  INTEGER :: jn, jgl, iodd, ik
854  INTEGER :: ins2, iallow, isym, iflag, itemax, jter
855 
856  REAL(KIND=JPRB) :: zhook_handle
857 
858  ! ------------------------------------------------------------------
859 
860  ! ------------------------------------------------------------------
861  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:LATITUDES_GAUSS',0,zhook_handle)
862  ! ------------------------------------------------------------------
863  !* 1. Initialization.
864  ! ---------------
865  !
866  !* 1.1 Calculation of ZFNLAT.
867  ! (Fourier coefficients of series expansion for
868  ! the ordinary Legendre polynomials).
869  ! Belousov, Swarztrauber use ZFN(0,0)=SQRT(2._JPRB)
870  ! IFS normalisation chosen to be 0.5*Integral(Pnm**2) = 1
871  zfn(0,0) = 2.
872  DO jn=1,kn
873  zfnn=zfn(0,0)
874  DO jgl=1,jn
875  zfnn = zfnn*sqrt(1.-0.25/(jgl**2))
876  ENDDO
877  iodd=mod(jn,2)
878  zfn(jn,jn)=zfnn
879  DO jgl=2,jn-iodd,2
880  zfn(jn,jn-jgl)=zfn(jn,jn-jgl+2)*float((jgl-1)*(2*jn-jgl+2))/float(jgl*(2*jn-jgl+1))
881  ENDDO
882  ENDDO
883 
884  iodd=mod(kn,2)
885  ik=iodd
886  DO jgl=iodd,kn,2
887  zfnlat(ik)=zfn(kn,jgl)
888  ik=ik+1
889  ENDDO
890 
891  !* 1.2 Find first approximation of the roots of the
892  ! Legendre polynomial of degree KN.
893  zpi=2.0*asin(1.0) ! Constant Pi
894  ins2 = kn/2+mod(kn,2)
895  DO jgl=1,ins2
896  z = (4*jgl-1)*zpi/(4*kn+2)
897  pl(jgl) = z+1.0/(tan(z)*8*(kn**2))
898  zreg(jgl) = cos(z)
899  zli(jgl) = cos(pl(jgl))
900  ENDDO
901 
902  ! ------------------------------------------------------------------
903 
904  !* 2. Computes roots and weights for transformed theta
905  ! ------------------------------------------------
906 
907  zeps = epsilon(z)
908  itemax = 20
909  iodd=mod(kn,2)
910 
911  DO jgl=ins2,1,-1
912 
913  ! * Initialization.
914 
915  zx = pl(jgl)
916  zdlx = zx
917  iflag = 0
918 
919  ! * Newton iteration.
920 
921  DO jter=1,itemax+1
922 
923  iter(jgl) = jter
924 
925  ! - Newton iteration step.
926 
927  zdlk = 0.0
928  IF( iodd==0 ) zdlk=0.5*zfnlat(0)
929  zzdlxn = 0.0
930  zdlldn = 0.0
931  ik=1
932 
933  IF(iflag == 0)THEN
934  DO jn=2-iodd,kn,2
935  ! normalised ordinary Legendre polynomial == \overbar{P_n}^0
936  zdlk = zdlk + zfnlat(ik)*cos(jn*zdlx)
937  ! normalised derivative == d/d\theta(\overbar{P_n}^0)
938  zdlldn = zdlldn - zfnlat(ik)*jn*sin(jn*zdlx)
939  ik=ik+1
940  ENDDO
941  ! Newton method
942  zdlmod = -zdlk/zdlldn
943  zzdlxn = zdlx+zdlmod
944  zxn = zzdlxn
945  zdlxn = zzdlxn
946  zmod(jgl) = zdlmod
947  ENDIF
948 
949  ! - Computes weight.
950 
951  IF(iflag == 1)THEN
952  DO jn=2-iodd,kn,2
953  ! normalised derivative
954  zdlldn = zdlldn - zfnlat(ik)*jn*sin(jn*zdlx)
955  ik=ik+1
956  ENDDO
957  zw = (2*kn+1)/zdlldn**2
958  ENDIF
959 
960  zx = zxn
961  zdlx = zdlxn
962 
963  IF(iflag == 1) EXIT
964  IF(abs(zmod(jgl)) <= zeps*1000.) iflag = 1
965  ENDDO
966 
967  pl(jgl) = zxn
968  pdl(jgl) = zdlxn
969  pw(jgl) = zw
970 
971  ENDDO
972 
973  ! convert to physical latitude space PMU
974  DO jgl=1,ins2
975  pl(jgl) = cos(pl(jgl))
976  pdl(jgl) = cos(pdl(jgl))
977  ENDDO
978 
979  DO jgl=1,kn/2
980  isym = kn-jgl+1
981  pl(isym) = -pl(jgl)
982  pdl(isym) = -pdl(jgl)
983  pw(isym) = pw(jgl)
984  ENDDO
985 
986  ! ------------------------------------------------------------------
987 
988  !* 3. Diagnostics.
989  ! ------------
990 
991  zra=6371229. ! Earth radius
992  DO jgl=1,ins2
993  zacos = acos(pl(jgl))
994  zm(jgl) = (zacos-acos(zli(jgl)))*zra
995  zr(jgl) = (zacos-acos(zreg(jgl)))*zra
996  zt(jgl) = zacos*180./zpi
997  ENDDO
998 
999  iallow = 10
1000  DO jgl=1,ins2
1001  IF(iter(jgl) > iallow)THEN
1002  WRITE(*,fmt='('' CONVERGENCE FAILED IN MODE_GRIDTYPE_GAUSS:LATITUDES_GAUSS '')')
1003  WRITE(*,fmt='('' ALLOWED : '',I4,''&
1004  &NECESSARY : '',&
1005  &I4)')iallow,iter(jgl)
1006  CALL abor1_sfx(' FAILURE IN MODE_GRIDTYPE_GAUSS:LATITUDES_GAUSS ')
1007  ENDIF
1008 
1009  !WRITE(*,FMT=&
1010  ! &'('' ROW ='',I4,'' ITERATIONS='',I4,'' ROOT='',F30.20,&
1011  ! &'' WEIGHT='',F30.20,'' MODIF :'',E8.2)')JGL,ITER(JGL),PL(JGL)&
1012  ! &,PW(JGL),PL(JGL)-ZLI(JGL)
1013  !WRITE(*,FMT=&
1014  ! &'(10X,'' LAST INC. : '',E8.2,'' MODIF IN M : '',F10.3,&
1015  ! &'' FROM THE REGULAR GRID : '',F10.3,'' COLAT '',F10.3)')&
1016  ! &ZMOD(JGL),ZM(JGL),ZR(JGL),ZT(JGL)
1017  ENDDO
1018 
1019  ! ------------------------------------------------------------------
1020 
1021  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:LATITUDES_GAUSS',1,zhook_handle)
1022  END SUBROUTINE latitudes_gauss
1023  !##################################
1024  !############################################################################
1025  !############################################################################
1026  !##################################
1027  SUBROUTINE mesh_size_gauss(KL,KNLATI,KNLOPA,PLAPO,PLOPO,PCODIL,&
1028  plat_xy,plat,plon,pmesh_size)
1029  !##################################
1030  !
1031  !!**** *MESH_SIZE_GAUSS* - routine to compute the global mesh size
1032  !!
1033  !! PURPOSE
1034  !! -------
1035  !!
1036  !!** METHOD
1037  !! ------
1038  !!
1039  !! AUTHOR
1040  !! ------
1041  !! *Meteo France*
1042  !!
1043  !! MODIFICATIONS
1044  !! -------------
1045  !! Original 06/2008
1046  !-------------------------------------------------------------------------------
1047  !
1048  USE modd_csts, ONLY : xradius, xpi
1049  !
1050  IMPLICIT NONE
1051  !
1052  INTEGER, INTENT(IN) :: kl ! number of point
1053  INTEGER, INTENT(IN) :: knlati ! number of pseudo-latitudes
1054  INTEGER, DIMENSION(KNLATI), INTENT(IN) :: knlopa ! number of pseudo-longitudes on each
1055  REAL, INTENT(IN) :: plapo ! latitude of the rotated pole (deg)
1056  REAL, INTENT(IN) :: plopo ! longitude of the rotated pole (deg)
1057  REAL, INTENT(IN) :: pcodil ! stretching factor (must be greater than or equal to 1)
1058  REAL, DIMENSION(KL), INTENT(IN) :: plat_xy ! pseudo latitudes
1059  REAL, DIMENSION(KL), INTENT(IN) :: plat ! latitude (degrees)
1060  REAL, DIMENSION(KL), INTENT(IN) :: plon ! longitude (degrees)
1061  REAL, DIMENSION(KL), INTENT(OUT):: pmesh_size ! global mesh size
1062  !
1063  INTEGER, DIMENSION(:), ALLOCATABLE :: ixx ! number of points in the latitude circle of each grid point
1064  INTEGER, DIMENSION(:), ALLOCATABLE :: iyy ! latitude circle of each grid point
1065  REAL, DIMENSION(:), ALLOCATABLE :: zdlat !
1066  REAL, DIMENSION(:), ALLOCATABLE :: plat_xy_c ! pseudo-latitude for each circle
1067  REAL, DIMENSION(:), ALLOCATABLE :: zxx ! X-length of each mesh in the gaussian grid
1068  REAL, DIMENSION(:), ALLOCATABLE :: zyy ! Y-length of each mesh in the gaussian grid
1069  REAL, DIMENSION(:), ALLOCATABLE :: zmap ! map factor
1070  INTEGER :: il
1071  INTEGER :: jl ! loop counter on number of points
1072  INTEGER :: jltot ! loop counter on total number of points
1073  REAL(KIND=JPRB) :: zhook_handle
1074  !
1075  !-------------------------------------------------------------------------------
1076  !
1077  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:MESH_SIZE_GAUSS',0,zhook_handle)
1078  ALLOCATE(ixx(kl))
1079  ALLOCATE(iyy(kl))
1080  ALLOCATE(zxx(kl))
1081  ALLOCATE(plat_xy_c(knlati))
1082  ALLOCATE(zdlat(knlati))
1083  ALLOCATE(zyy(kl))
1084  ALLOCATE(zmap(kl))
1085  !
1086  ! 1.0 Length according to the X axis (longitudes)
1087  ! ------------------------------
1088  !
1089  il=1
1090  DO jl = 1,knlati
1091  DO jltot = 1,knlopa(jl)
1092  ixx(il)=knlopa(jl)
1093  iyy(il)=jl
1094  il=il+1
1095  ENDDO
1096  ENDDO
1097  !
1098  !
1099  DO jl = 1,kl
1100  zxx(jl) = 2.*xpi*xradius*cos(plat_xy(jl)*xpi/180.)/float(ixx(jl))
1101  ENDDO
1102  !
1103  ! 2.0 Length according to the Y axis (latitudes)
1104  ! ------------------------------
1105  !
1106  il=1
1107  DO jl = 1,knlati
1108  plat_xy_c(jl)=plat_xy(il)
1109  il=il+knlopa(jl)
1110  ENDDO
1111  !
1112  zdlat(1)=90.-plat_xy_c(1)+((plat_xy_c(1)-plat_xy_c(2))/2.)
1113  DO jl = 2,knlati-1
1114  zdlat(jl)=((plat_xy_c(jl-1)-plat_xy_c(jl))/2.)+((plat_xy_c(jl)-plat_xy_c(jl+1))/2.)
1115  ENDDO
1116  zdlat(knlati)=((plat_xy_c(knlati-1)-plat_xy_c(knlati))/2.)+plat_xy_c(knlati)+90.
1117  zdlat(:)=zdlat(:)*xpi*xradius/180.
1118  !
1119  DO jl = 1,kl
1120  zyy(jl)=zdlat(iyy(jl))
1121  ENDDO
1122  !
1123  ! 3.0 grid mesh
1124  ! ---------
1125  !
1126  CALL map_factor_gauss(plapo,plopo,pcodil,plat,plon,zmap)
1127  !
1128  pmesh_size(:) = zxx(:) * zyy(:) * zmap(:)**2
1129  !
1130  DEALLOCATE(zxx )
1131  DEALLOCATE(zyy )
1132  DEALLOCATE(zmap)
1133  DEALLOCATE(plat_xy_c)
1134  DEALLOCATE(zdlat)
1135  DEALLOCATE(ixx )
1136  DEALLOCATE(iyy )
1137  IF (lhook) CALL dr_hook('MODE_GRIDTYPE_GAUSS:MESH_SIZE_GAUSS',1,zhook_handle)
1138  !
1139  !-------------------------------------------------------------------------------
1140  END SUBROUTINE mesh_size_gauss
1141 
1142  !############################################################################
1143  !############################################################################
1144  !############################################################################
1145 
1146 END MODULE mode_gridtype_gauss
1147 
mode_gridtype_gauss
Definition: mode_gridtype_gauss.F90:6
mode_geo_gauss::map_fac
Definition: mode_geo_gauss.F90:48
mode_gridtype_gauss::put_gridtype_gauss
subroutine put_gridtype_gauss(PGRID_PAR, KNLATI, PLAPO, PLOPO, PCODIL, KNLOPA, KL, PLAT, PLON, PLAT_XY, PLON_XY, PMESH_SIZE, PLONINF, PLATINF, PLONSUP, PLATSUP)
Definition: mode_gridtype_gauss.F90:26
mode_gridtype_gauss::get_gridtype_gauss
subroutine get_gridtype_gauss(PGRID_PAR, KNLATI, PLAPO, PLOPO, PCODIL, KNLOPA, KL, PLAT, PLON, PLAT_XY, PLON_XY, PMESH_SIZE, PLONINF, PLATINF, PLONSUP, PLATSUP)
Definition: mode_gridtype_gauss.F90:107
mode_gridtype_gauss::gauss_grid_limits
subroutine gauss_grid_limits(KNLATI, KNLOPA, PXINF, PXSUP, PYINF, PYSUP)
Definition: mode_gridtype_gauss.F90:458
eggangles::angle_domain
Definition: dummy_eggangles.F90:11
abor1_sfx
subroutine abor1_sfx(YTEXT)
Definition: abor1_sfx.F90:6
eggangles::p_asin
Definition: dummy_eggangles.F90:35
mode_geo_gauss
Definition: mode_geo_gauss.F90:5
mode_gridtype_gauss::mesh_size_gauss
subroutine mesh_size_gauss(KL, KNLATI, KNLOPA, PLAPO, PLOPO, PCODIL, PLAT_XY, PLAT, PLON, PMESH_SIZE)
Definition: mode_gridtype_gauss.F90:1027
mode_gridtype_gauss::latlon_gauss
subroutine latlon_gauss(PLON_XY, PLAT_XY, KL, PLOPO, PLAPO, PCODIL, PLON, PLAT)
Definition: mode_gridtype_gauss.F90:252
mode_gridtype_gauss::xy_gauss
subroutine xy_gauss(PCODIL, KSIZE, PNODATA, PVALUE, PLAT_XY, PLON_XY)
Definition: mode_gridtype_gauss.F90:551
mode_gridtype_gauss::latitudes_gauss
subroutine latitudes_gauss(KN, PL, PDL, PW)
Definition: mode_gridtype_gauss.F90:764
modd_get_mesh_index_gauss
Definition: modd_get_mesh_index_gauss.F90:6
mode_gridtype_gauss::map_factor_gauss
subroutine map_factor_gauss(PLAPO, PLOPO, PCODIL, PLAT, PLON, PMAP)
Definition: mode_gridtype_gauss.F90:683
mode_gridtype_gauss::comp_gridtype_gauss
subroutine comp_gridtype_gauss(KNLATI, KNLOPA, KL, KTYP, PLAT_XY, PLON_XY)
Definition: mode_gridtype_gauss.F90:367
modd_csts
Definition: modd_csts.F90:6
Generated on Fri May 13 2016 09:31:04 for SURFEX by   doxygen 1.8.5