source: LMDZ.3.3/trunk/libf/dyn3d/fxyhypF @ 507

Last change on this file since 507 was 2, checked in by lmdz, 25 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 16.1 KB
Line 
1       SUBROUTINE fxyhyp( clat,  rlatu,yprimu ,  rlatv,yprimv   ,
2     *                         rlatu1 ,yprimu1, rlatu2,yprimu2   )
3c
4       IMPLICIT NONE
5c   -----------------------
6c     Auteur  :  P. Le Van
7c   -----------------------
8c
9c .. Objet : Calcul d'une fonction f(y)  ayant pour derivee une tangente
10c        hyperbolique , donc des ordonnees y de la grille horizontale ,
11c       ainsi que ses derivees yprim , utlisees dans  Inigeom  .
12c
13c
14c ....  rlatu   =  fxyhyp (   j      )     avec 1  <= j  <= jjm + 1
15c ....  rlatv   =  fxyhyp ( j + 0.5  )     avec 1  <= j  <= jjm
16c ....  rlatu1  =  fxyhyp ( j + 0.25 )     avec 1  <= j  <= jjm
17c ....  rlatu2  =  fxyhyp ( j + 0.75 )     avec 1  <= j  <= jjm
18c
19
20       INTEGER nmax
21       PARAMETER ( nmax = 100 000 )
22c
23#include "dimensions.h"
24#include "paramet.h"
25#include "comconst.h"
26c
27c    .................   Arguments     ................................
28c
29c   clat est (en degres) la latitude du centre du zoom : arg. d'entree
30c    Les autres arguments     sont                des argum. de sortie
31c
32       REAL clat
33       REAL rlatu(jjp1),rlatv(jjm),yprimu(jjp1),yprimv(jjm)
34       REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm)
35c
36c    .................................................................
37c
38c             ...... .  Variables  locales   .........   
39c
40       REAL yprim( 0:nmax ),yi( 0:nmax ),yf( 0:nmax )
41       REAL yyy(0:nmax)
42       REAL delta,den,eps,phij,pijjm,pis2,pis2pyo,pism,psi,psi0
43       REAL rappis2,unpdel,y,yint,yj,yo,yo1,yppp,ypr
44       REAL yprimin,pisnmax,ydeg
45       INTEGER i,j,it,iter
46c
47c
48c -------------------------------------------------------------------------
49   
50c   On veut calculer  y(Y) et y'(Y) , avec  Y ,latitude de travail =f(j),
51c       et   y , latitude vraie .
52c
53c  Pour cela ,on se donne pour  Y'(y) une fonction hyperbolique calculee
54c  sur  ( 1+ nmax ) points ( plus nmax est grand, plus la precision pour
55c  y(Y)  sera grande ) , on integre Y'(y) sur ces ( 1+nmax) pts  pour
56c  obtenir  Y(y) apres normalisation  , puis on  aura y(Y)  en resolvant
57c   l'equation  Y(yj) = Yj  par Newton_Raphson ,  avec Yj = f(j) =
58c   pi/2 + pi/jjm ( 1.- j)  pour les scalaires  ( j =1, jjm+1 )
59c   On en tire donc  yj  qui correspond  a  la ligne  j   .
60c   On a ensuite pour les derivees :  y'(Y) = 1./ Y'(y) * ( pi/ jjm )
61c
62c   On calculera   y(Y) et y'(Y)  successivement   aux  latitudes  de  U
63c     ,  V  ,  U + 0.25 , U + 0.75   .
64c
65c   Notations ;
66c  -------------
67c     yi    represente les nmax+1 latitudes de depart ( -pi/2 a pi/2 )
68c  yyy(nmax)  repres.  l'integrale de Y'(y) de - pi/2 a pi/2
69c    yprim  represente  Y'(y)
70c     yf    represente  Y (y)
71c
72c     yo    represente la latitude du centre du zoom
73c    delta  represente la demi largeur du zoom
74c
75c
76c  N.B :  La fonction hyperbolique Y'(y) est  definie comme suit :
77c ------
78c
79c    Y' = TANH ( ( psio - psi )/psi*(1- psi) )  avec  :
80c
81c          psi =  ( y - yo ) / ( pi/2 - yo )  si     yo < y < pi/2
82c          psi =  (yo -  y )/  ( pi/2 + yo )  si - pi/2 < y < yo
83c
84c
85c  ....  yo    represente la latitude du centre du zoom         .....
86c
87c  ....  psi0  et delta sont des parametres variant entre 0 et 1  ....
88c----------------------------------------------------------------------
89c
90cc     psi0 = 0.3 et  delta = 0.5  sont les valeurs qui conviennent
91c      pour  *** iim = 96  et jjm = 72 ***  , d'apres des tests sur
92c      les valeurs du Laplacien itere ( 2 fois) , analytiques et cal-
93c      culees.   Pour d'autres valeurs de  iim et jjm , on peut tester
94c      les nouvelles valeurs de psi0 et delta en lancant  testfxyhyp .
95c
96c ----------------------------------------------------------------------
97c    Fxyhp  est appele par inigeom  a la place de  fxynew (  fonction
98c    a derivee sinusoidale )  si la variable logique fxyhypb , lue par
99c    defrun_new  est  . TRUE.    .
100c ----------------------------------------------------------------------
101c
102
103       pi     =  2.*ASIN( 1. )
104       pis2   =  pi/2.
105       yo     =  clat  * pi/180.
106       psi0   =  0.3
107       delta  =  0.5
108       unpdel =  1. + delta
109       eps    =  .1e-7
110c
111       PRINT *,' ----------------------------------------------------'
112       PRINT 1,    psi0,delta,clat
113       PRINT *,' ----------------------------------------------------'
114
115       pisnmax =  pi/FLOAT(nmax)
116       yi(0)   = - pis2
117
118       DO i = 1, nmax
119        yi(i)= - pis2 + FLOAT(i)*pisnmax
120       ENDDO
121
122       yyy(0)=  0.
123
124c
125c     ***************************************************************
126c     ************      Cas   ou   yo  <  0.     ********************
127c     ***************************************************************
128c
129       IF( yo.LT. 0. )   THEN
130c
131       pis2pyo = - pis2 + yo
132       rappis2 = ( -pis2-yo)/pis2pyo
133       rappis2 = rappis2 * rappis2
134c
135        DO 5 i = 1, nmax
136
137        y = 0.5 * ( yi(i)+ yi(i-1) )
138        IF( y.LT.- pis2.OR.y.GT.pis2 )  STOP 0
139
140        IF(y.GT.yo)  THEN
141        psi = (y-yo)/(pis2-yo)
142        ypr =  TANH( (psi0 -psi)/(psi*(1.-psi)) )
143         
144        ELSE
145        psi = (yo-y)/(pis2+yo)
146        ypr = 1.-rappis2 + rappis2* TANH ( (psi0 -psi)/(psi*(1.-psi)) )
147        ENDIF
148
149        yyy(i)= yyy(i-1)+ ypr*(yi(i)-yi(i-1))
150
151 5    CONTINUE
152
153       den = unpdel + yyy(nmax)/pi
154
155       yyy(0)=   0.
156       yf(0) = - pis2
157
158       yprim(nmax) = -1.
159       yprim(0)    = 1.- rappis2 - rappis2
160
161       DO 12 i =  0 ,nmax
162        y = yi(i)
163       IF(y.LT.-pis2.OR.y.GT.pis2)  STOP 1
164
165       IF(y.LT.yo)  THEN
166           psi = (yo -y)/(pis2+yo)
167         IF(i.NE.0.AND.i.NE.nmax)  THEN
168         yprim(i) = 1.- rappis2 + rappis2
169     *                 * TANH( (psi0 -psi)/(psi*(1.-psi)) )
170         ENDIF
171       ELSE
172        psi = (y-yo)/(pis2-yo)
173          IF(i.NE.0.AND.i.NE.nmax) THEN
174           yprim(i) =  TANH( (psi0 -psi)/(psi*(1.-psi)) )
175          ENDIF
176       ENDIF
177
178c
179c   .....................................................
180c   .......   Normalisation de  Y'(y) et de  Y(y)   .....
181c   .....................................................
182c
183        yprim(i)= ( yprim(i)+ unpdel )/ den
184        yf(i)= (yyy(i)+  unpdel*( yi(i)+ pis2))/ den  - pis2
185c
186 12    CONTINUE
187c
188      ELSE
189c
190c     **************************************************************
191c     **********          Cas  ou  yo  >=  0 .        **************
192c     **************************************************************
193c
194       pis2pyo = pis2 + yo
195       rappis2 = ( pis2 - yo )/pis2pyo
196       rappis2 = rappis2 * rappis2
197c
198
199       DO 13 i = 1, nmax
200        y= 0.5 * ( yi(i)+ yi(i-1) )
201        IF(y.LT.- pis2.OR.y.GT.pis2)   STOP 2
202
203       IF(y.LT.yo)  THEN
204        psi = (yo -y)/(pis2+yo)
205        ypr =  TANH( (psi0 - psi)/(psi*(1.-psi)) )
206
207       ELSE
208        psi = (y- yo)/(pis2-yo)
209        ypr = 1.-rappis2 + rappis2* TANH ( (psi0 -psi)/(psi*(1.-psi)) )
210
211       ENDIF
212
213        yyy(i)= yyy(i-1)+ ypr*(yi(i)-yi(i-1))
214
215 13    CONTINUE
216
217       den = unpdel + yyy(nmax)/pi
218
219       yyy(0)= 0.
220       yf(0) = -pis2
221
222       yprim(0)   = -1.
223       yprim(nmax)= 1.- rappis2 - rappis2
224
225       DO 15 i = 0,nmax
226        y = yi(i)
227       IF(y.LT.-pis2.OR.y.GT.pis2)   STOP 3
228
229       IF(y.LT.yo)  THEN
230        psi = (yo -y)/(pis2+yo)
231        IF(i.NE.0.AND.i.NE.nmax)
232     *   yprim(i) =  TANH( (psi0 -psi)/(psi*(1.-psi)) )
233       ELSE
234        psi = (y- yo)/(pis2-yo)
235        IF(i.NE.0.AND.i.NE.nmax)
236     *   yprim(i) = 1.- rappis2 + rappis2
237     *               * TANH( (psi0 -psi)/(psi*(1.-psi)) )                             
238       ENDIF
239c   .....................................................
240c
241c   .......   Normalisation de  Y'(y) et de  Y(y)   .....
242c   .....................................................
243c
244        yprim(i)= ( yprim(i)+ unpdel )/ den
245        yf(i)= (yyy(i)+  unpdel*( yi(i)+ pis2))/ den  - pis2
246
247 15    CONTINUE
248
249       ENDIF
250c
251c     ***************************************************************
252
253      yf(0)   = - pis2
254      yf(nmax)=   pis2
255       
256      pijjm = pi/FLOAT(jjm)
257
258c   ...............................................................
259c
260c     Calculs  de y(Y) et y'(Y)   pour les latitudes  de  U   .....
261c
262      PRINT *,' ***************************** '
263      PRINT *,' ***  Calcul de rlatu    ***** '
264      PRINT *,' ***************************** '
265c
266c   ...............................................................
267c
268      PRINT *,' *** j   rlatu  yprimu **** '
269c
270      DO 1500 j=1, jjm +1
271
272      phij = pis2 + pijjm* ( 1.-Float(j) )
273      yo1  = 0.
274      yj   = yo1
275c
276      DO 500 iter = 1,300
277
278      DO 250 it = nmax,0,-1
279      IF( yj.GE.yi(it))  GO TO 350
280250   CONTINUE
281
282      it= 0
283      yj = yi(it)
284
285350   CONTINUE
286       IF(it.EQ.nmax)  THEN
287         yint    = yf(nmax)
288         yprimin = yprim(nmax)
289      ELSE
290c  .................................................................
291c  ....  Interpolation entre  yi(it) et yi(it+1)   pour avoir Y(yj) 
292c      .....           et   Y'(yj)                             .....
293c  .................................................................
294      yint = (yf(it+1)-yf(it))/(yi(it+1)-yi(it))* ( yj-yi(it) )     +
295     +                           yf(it)
296      yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) +
297     +                           yprim(it)
298      ENDIF
299      yj = yj - (yint-phij)/yprimin
300
301      IF( ABS(yj-yo1).LE.eps) GO TO 550
302      yo1 = yj
303c
304500   CONTINUE
305      PRINT *,' ***   PAS DE SOLUTION  ****  ',j,phij
306        STOP 4
307550   CONTINUE
308       ydeg     = yj * 180./pi
309       rlatu(j) = yj
310       IF(j.EQ.1) rlatu(j) = pis2
311c
312c    .....  calcul de  yprimu .....
313c
314
315       IF(j.EQ.1.OR.j.EQ.jjm +1) GO TO 1500
316       IF(yj.LT.yo)  THEN
317        psi = (yo -yj)/(pis2+yo)
318       ELSE
319        psi = (yj-yo)/(pis2-yo)
320       ENDIF
321
322      DO  it = nmax,0,-1
323      IF( yj.GE.yi(it))  GO TO 555
324      ENDDO
325555    CONTINUE
326       IF(it.EQ.nmax) THEN
327       PRINT *,' IT = nmax  Pbs ..',it,j
328       STOP 5
329       ENDIF
330       yppp=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it)) +
331     *               yprim(it)
332       yprimu(j) = pijjm/yppp
333c
334       PRINT 3000,j,ydeg,yprimu(j)
3351500   CONTINUE
336       rlatu(1)     =   pis2
337       rlatu(jjm+1) = - pis2
338c
339c     .............................................
340c
341c     .....   Calculs  pour  rlatv et yprimv  .....
342c     .............................................
343c
344
345      PRINT *,' ***************************** '
346      PRINT *,' ****  Calcul  de rlatv  ***** '
347      PRINT *,' ***************************** '
348      PRINT *,' *** j   rlatv  yprimv ****  '
349
350      DO 1600 j=1, jjm
351
352      phij  = pis2+ pijjm*(1.-(FLOAT(j)+0.5))
353      yo1   = 0.
354      yj    = yo1
355
356      DO 1550 iter =1,300
357
358      DO 1520 it = nmax,0,-1
359      IF( yj.GE.yi(it))  GO TO 1530
3601520   CONTINUE
361c
362      it  = 0
363      yj  = yi(it)
3641530   CONTINUE
365       IF(it.EQ.nmax)  THEN
366       yint    = yf(nmax)
367       yprimin = yprim(nmax)
368      ELSE
369      yint= (yf(it+1)-yf(it))/(yi(it+1)-yi(it)) * ( yj-yi(it) )  +
370     +                          yf(it)
371      yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) +
372     +                          yprim(it)
373      ENDIF
374      yj = yj - (yint-phij)/yprimin
375
376      IF( ABS(yj-yo1).LE.eps ) GO TO 1560
377      yo1 = yj
378c
3791550   CONTINUE
380      PRINT *,'  ***  PAS DE SOLUTION  ***  ',j,phij
3811560   CONTINUE
382      ydeg    = yj * 180./pi
383      rlatv(j)= yj
384c
385c    .....  calcul de  yprimv .....
386c
387       IF(yj.LT.-pis2.OR.yj.GT.pis2)   STOP 6
388
389       IF(yj.LT.yo)  THEN
390        psi = (yo -yj)/(pis2+yo)
391       ELSE
392        psi = (yj-yo)/(pis2-yo)
393       ENDIF
394      DO  it = nmax,0,-1
395      IF( yj.GE.yi(it))   GO TO 1570
396      ENDDO
3971570    CONTINUE
398       IF(it.EQ.nmax) THEN
399       PRINT *,' IT = nmax  Pbs ..',it,j
400       STOP 7
401       ENDIF
402       yppp = (yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it)) +
403     *                   yprim(it)
404       yprimv(j)= pijjm/ yppp
405c
406       PRINT 3000,j,ydeg,yprimv(j)
407
4081600  CONTINUE
409
410c
411c    ..................................................
412c
413      PRINT *,' ************************************ '
414      PRINT *,' ****  Calcul pour f( j + 0.25 ) **** '
415      PRINT *,' ************************************ '
416      PRINT *,' *** j  rlatu1  yprimu1 **** '
417c
418c    ..................................................
419c
420      DO 1800 j=1, jjm
421
422      phij = pis2+ pijjm*( 1.- (FLOAT(j) + 0.25) )
423      yo1  = 0.
424      yj   = yo1
425C
426
427      DO 1620 iter = 1, 500
428
429      DO 1610 it = nmax,0,-1
430      IF( yj.GE.yi(it))  GO TO 1615
4311610   CONTINUE
432c
433      it = 0
434      yj = yi(it)
4351615   CONTINUE
436      IF(it.EQ.nmax)  THEN
437       yint    = yf(nmax)
438       yprimin = yprim(nmax)
439      ELSE
440      yint= (yf(it+1)-yf(it))/(yi(it+1)-yi(it))*(yj-yi(it))+
441     * yf(it)
442      yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it))+
443     * yprim(it)
444      ENDIF
445      yj = yj - (yint-phij)/yprimin
446c      IF( (- pis2-yj)*(yj- pis2 ).LT.0.)PRINT *,' Newton hors limites'
447
448      IF( ABS(yj-yo1).LE.eps ) GO TO 1625
449      yo1 = yj
450c
4511620   CONTINUE
452      PRINT *,' ***  PAS DE SOLUTION  ****  ',j,phij
4531625   CONTINUE
454      ydeg     = yj * 180./pi
455      rlatu1(j)= yj
456c
457c    .....  calcul de  yprimu1 .....
458c
459       IF( yj.LT.-pis2.OR.yj.GT.pis2 )  STOP 8
460c
461       IF(yj.LT.yo)  THEN
462        psi = (yo -yj)/(pis2+yo)
463       ELSE
464        psi = (yj-yo)/(pis2-yo)
465       ENDIF
466      DO  it = nmax,0,-1
467      IF( yj.GE.yi(it))   GO TO 1630
468      ENDDO
4691630    CONTINUE
470       IF(it.EQ.nmax) THEN
471       PRINT *,' IT = nmax  Pbs ..',it,j
472       STOP 9
473       ENDIF
474       yppp = (yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it)) +
475     *          yprim(it)
476       yprimu1(j) = pijjm/yppp
477c
478       PRINT 3000,j,ydeg,yprimu1(j)
479
4801800  CONTINUE
481   
482c    ..............................................
483c
484      PRINT *,' ********************************** '
485      PRINT *,' **** Calcul pour f( j + 3/4  ) *** '
486      PRINT *,' ********************************** '
487      PRINT *,' *** j  rlatu2  yprimu2 ***  '
488c    ..............................................
489c
490      DO 1900 j=1, jjm
491
492      phij = pis2+ pijjm*(1.-(FLOAT(j)+0.75))
493      yo1  = 0.
494      yj   = yo1
495
496      DO 1700 iter =1,500
497
498      DO 1680 it = nmax,0,-1
499      IF( yj.GE.yi(it))   GO TO 1690
5001680   CONTINUE
501      it  = 0
502      yj  = yi(it)
5031690   CONTINUE
504      IF(it.EQ.nmax)  THEN
505       yint     = yf(nmax)
506       yprimin = yprim(nmax)
507      ELSE
508      yint= (yf(it+1)-yf(it))/(yi(it+1)-yi(it)) * ( yj -yi(it) ) +
509     +                             yf(it)
510      yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) +
511     +                             yprim(it)
512      ENDIF
513      yj = yj - (yint-phij)/yprimin
514
515      IF( ABS(yj-yo1).LE.eps ) GO TO 1705
516      yo1 = yj
517c
5181700   CONTINUE
519      PRINT *,' PAS DE SOLUTION ',j,phij
520        STOP 10
5211705   CONTINUE
522      ydeg      = yj * 180./pi
523      rlatu2(j) = yj
524c
525c    .....  calcul de  yprimu2 .....
526c
527       IF(yj.LT.- pis2.OR.yj.GT.pis2)  STOP 11
528
529       IF(yj.LT.yo)  THEN
530        psi = (yo -yj)/(pis2+yo)
531       ELSE
532        psi = (yj-yo)/(pis2-yo)
533       ENDIF
534c
535       DO  it = nmax,0,-1
536       IF( yj.GE.yi(it))  GO TO 1710
537       ENDDO
5381710    CONTINUE
539       IF(it.EQ.nmax) THEN
540       PRINT *,' IT = nmax  Pbs ..',it,j
541        STOP 12
542       ENDIF
543      yppp =(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj - yi(it) ) +
544     +                       yprim(it)
545       yprimu2(j) = pijjm/yppp
546c
547       PRINT 3000,j,ydeg,yprimu2(j)
548
5491900  CONTINUE
550c
551       PRINT *,'  *** TEST DE COHERENCE DANS  FXYP **** '
552c
553       DO j = 1, jjm
554c
555       IF(rlatu1(j).LE.rlatu2(j))   THEN
556         PRINT *,' Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j
557         STOP 13
558       ENDIF
559c
560       IF(rlatu2(j).LE.rlatu(j+1))  THEN
561        PRINT *,' Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j
562        STOP 14
563       ENDIF
564c
565       IF(rlatu(j).LE.rlatu1(j))    THEN
566        PRINT *,' rlatu < rlatu1 ',rlatu(j),rlatu1(j),j
567        STOP 15
568       ENDIF
569c
570       IF(rlatv(j).LE.rlatu2(j))    THEN
571        PRINT *,' rlatv > rlatu2 ',rlatv(j),rlatu2(j),j
572        STOP 16
573       ENDIF
574c
575       IF(rlatv(j).ge.rlatu1(j))    THEN
576        PRINT *,' rlatv < rlatu1 ',rlatv(j),rlatu1(j),j
577        STOP 17
578       ENDIF
579c
580       IF(rlatv(j).ge.rlatu(j))     THEN
581        PRINT *,' rlatv > rlatu ',rlatv(j),rlatu(j),j
582        STOP 18
583       ENDIF
584c
585       ENDDO
586c
587       PRINT *,' ***  Les latitudes yprimu,yprimv,yprimu1,yprimu2 ',
588     * ' sont coherentes entre elles !  *** '
589
590c
591 1     FORMAT( 3x,'Fxyhyp.  psi0  delta clat ',3f7.2)
5923000   FORMAT(2x,i4,2x,f8.2,f8.3)
593       RETURN
594       END
Note: See TracBrowser for help on using the repository browser.