source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn3d_common/fyhyp.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 11.6 KB
Line 
1!
2! $Id: fyhyp.F 1403 2010-07-01 09:02:53Z fairhead $
3!
4c
5c
6       SUBROUTINE fyhyp ( yzoomdeg, grossism, dzooma,tau  , 
7     ,  rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
8     ,  champmin,champmax                                            ) 
9
10cc    ...  Version du 01/04/2001 ....
11
12       IMPLICIT NONE
13c
14c    ...   Auteur :  P. Le Van  ... 
15c
16c    .......    d'apres  formulations  de R. Sadourny  .......
17c
18c     Calcule les latitudes et derivees dans la grille du GCM pour une
19c     fonction f(y) a tangente  hyperbolique  .
20c
21c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc)
22c     dzoom  etant  la distance totale de la zone du zoom ( en radians )
23c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom   
24c
25c
26c N.B : Il vaut mieux avoir : grossism * dzoom  <  pi/2  (radians) ,en lati.
27c      ********************************************************************
28c
29c
30!-----------------------------------------------------------------------
31!   INCLUDE 'dimensions.h'
32!
33!   dimensions.h contient les dimensions du modele
34!   ndm est tel que iim=2**ndm
35!-----------------------------------------------------------------------
36
37      INTEGER iim,jjm,llm,ndm
38
39      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
40
41!-----------------------------------------------------------------------
42!
43! $Header$
44!
45!
46!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
47!                 veillez  n'utiliser que des ! pour les commentaires
48!                 et  bien positionner les & des lignes de continuation
49!                 (les placer en colonne 6 et en colonne 73)
50!
51!
52!-----------------------------------------------------------------------
53!   INCLUDE 'paramet.h'
54
55      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
56      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
57      INTEGER  ijmllm,mvar
58      INTEGER jcfil,jcfllm
59
60      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
61     &    ,jjp1=jjm+1-1/jjm)
62      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
63      PARAMETER( kftd  = iim/2 -ndm )
64      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
65      PARAMETER( ip1jmi1= ip1jm - iip1 )
66      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
67      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
68      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
69
70!-----------------------------------------------------------------------
71
72       INTEGER      nmax , nmax2
73       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
74c
75c
76c     .......  arguments  d'entree    .......
77c
78       REAL yzoomdeg, grossism,dzooma,tau
79c         ( rentres  par  run.def )
80
81c     .......  arguments  de sortie   .......
82c
83       REAL rrlatu(jjp1), yyprimu(jjp1),rrlatv(jjm), yyprimv(jjm),
84     , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
85
86c
87c     .....     champs  locaux    .....
88c
89     
90       REAL   dzoom
91       REAL(KIND=8) ylat(jjp1), yprim(jjp1)
92       REAL(KIND=8) yuv
93       REAL(KIND=8) yt(0:nmax2)
94       REAL(KIND=8) fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2)
95       SAVE Ytprim, yt,Yf
96       REAL(KIND=8) Yf(0:nmax2),yypr(0:nmax2)
97       REAL(KIND=8) yvrai(jjp1), yprimm(jjp1),ylatt(jjp1)
98       REAL(KIND=8) pi,depi,pis2,epsilon,y0,pisjm
99       REAL(KIND=8) yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax
100       REAL(KIND=8) yfi,Yf1,ffdy
101       REAL(KIND=8) ypn,deply,y00
102       SAVE y00, deply
103
104       INTEGER i,j,it,ik,iter,jlat
105       INTEGER jpn,jjpn
106       SAVE jpn
107       REAL(KIND=8) a0,a1,a2,a3,yi2,heavyy0,heavyy0m
108       REAL(KIND=8) fa(0:nmax2),fb(0:nmax2)
109       REAL y0min,y0max
110
111       REAL(KIND=8)     heavyside
112
113       pi       = 2. * ASIN(1.)
114       depi     = 2. * pi
115       pis2     = pi/2.
116       pisjm    = pi/ REAL(jjm)
117       epsilon  = 1.e-3
118       y0       =  yzoomdeg * pi/180.
119
120       IF( dzooma.LT.1.)  THEN
121         dzoom = dzooma * pi
122       ELSEIF( dzooma.LT. 12. ) THEN
123         WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug
124     ,menter et relancer ! '
125         STOP 1
126       ELSE
127         dzoom = dzooma * pi/180.
128       ENDIF
129
130       WRITE(6,18)
131       WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)'
132       WRITE(6,24) y0,grossism,tau,dzoom
133
134       DO i = 0, nmax2 
135        yt(i) = - pis2  + REAL(i)* pi /nmax2
136       ENDDO
137
138       heavyy0m = heavyside( -y0 )
139       heavyy0  = heavyside(  y0 )
140       y0min    = 2.*y0*heavyy0m - pis2
141       y0max    = 2.*y0*heavyy0  + pis2
142
143       fa = 999.999
144       fb = 999.999
145       
146       DO i = 0, nmax2 
147        IF( yt(i).LT.y0 )  THEN
148         fa (i) = tau*  (yt(i)-y0+dzoom/2. )
149         fb(i) =   (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) )
150        ELSEIF ( yt(i).GT.y0 )  THEN
151         fa(i) =   tau *(y0-yt(i)+dzoom/2. )
152         fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 ) 
153       ENDIF
154       
155       IF( 200.* fb(i) .LT. - fa(i) )   THEN
156         fhyp ( i) = - 1.
157       ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
158         fhyp ( i) =   1.
159       ELSE 
160         fhyp(i) =  TANH ( fa(i)/fb(i) )
161       ENDIF
162
163       IF( yt(i).EQ.y0 )  fhyp(i) = 1.
164       IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1.
165
166       ENDDO
167
168cc  ....  Calcul  de  beta  ....
169c
170       ffdy   = 0.
171
172       DO i = 1, nmax2
173        ymoy    = 0.5 * ( yt(i-1) + yt( i ) )
174        IF( ymoy.LT.y0 )  THEN
175         fa(i)= tau * ( ymoy-y0+dzoom/2.) 
176         fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy )
177        ELSEIF ( ymoy.GT.y0 )  THEN
178         fa(i)= tau * ( y0-ymoy+dzoom/2. ) 
179         fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 )
180        ENDIF
181
182        IF( 200.* fb(i) .LT. - fa(i) )    THEN
183         fxm ( i) = - 1.
184        ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
185         fxm ( i) =   1.
186        ELSE
187         fxm(i) =  TANH ( fa(i)/fb(i) )
188        ENDIF
189         IF( ymoy.EQ.y0 )  fxm(i) = 1.
190         IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1.
191         ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) )
192
193        ENDDO
194
195        beta  = ( grossism * ffdy - pi ) / ( ffdy - pi )
196
197       IF( 2.*beta - grossism.LE. 0.)  THEN
198
199        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
200     ,tine fyhyp est mauvaise ! '
201        WRITE(6,*)'Modifier les valeurs de  grossismy ,tauy ou dzoomy',
202     , ' et relancer ! ***  '
203        CALL ABORT_GCM("FYHYP", "", 1)
204
205       ENDIF
206c
207c   .....  calcul  de  Ytprim   .....
208c
209       
210       DO i = 0, nmax2
211        Ytprim(i) = beta  + ( grossism - beta ) * fhyp(i)
212       ENDDO
213
214c   .....  Calcul  de  Yf     ........
215
216       Yf(0) = - pis2
217       DO i = 1, nmax2
218        yypr(i)    = beta + ( grossism - beta ) * fxm(i)
219       ENDDO
220
221       DO i=1,nmax2
222        Yf(i)   = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) )
223       ENDDO
224
225c    ****************************************************************
226c
227c   .....   yuv  = 0.   si calcul des latitudes  aux pts.  U  .....
228c   .....   yuv  = 0.5  si calcul des latitudes  aux pts.  V  .....
229c
230      WRITE(6,18)
231c
232      DO 5000  ik = 1,4
233
234       IF( ik.EQ.1 )  THEN
235         yuv  = 0.
236         jlat = jjm + 1
237       ELSE IF ( ik.EQ.2 )  THEN
238         yuv  = 0.5
239         jlat = jjm 
240       ELSE IF ( ik.EQ.3 )  THEN
241         yuv  = 0.25
242         jlat = jjm 
243       ELSE IF ( ik.EQ.4 )  THEN
244         yuv  = 0.75
245         jlat = jjm 
246       ENDIF
247c
248       yo1   = 0.
249       DO 1500 j =  1,jlat
250        yo1   = 0.
251        ylon2 =  - pis2 + pisjm * ( REAL(j)  + yuv  -1.) 
252        yfi    = ylon2
253c
254       DO 250 it =  nmax2,0,-1
255        IF( yfi.GE.Yf(it))  GO TO 350
256250    CONTINUE
257       it = 0
258350    CONTINUE
259
260       yi = yt(it)
261       IF(it.EQ.nmax2)  THEN
262        it       = nmax2 -1
263        Yf(it+1) = pis2
264       ENDIF
265c  .................................................................
266c  ....  Interpolation entre  yi(it) et yi(it+1)   pour avoir Y(yi) 
267c      .....           et   Y'(yi)                             .....
268c  .................................................................
269
270       CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1),   
271     ,                  yt(it),yt(it+1) ,   a0,a1,a2,a3   )     
272
273       Yf1     = Yf(it)
274       Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi
275
276       DO 500 iter = 1,300
277         yi = yi - ( Yf1 - yfi )/ Yprimin
278
279        IF( ABS(yi-yo1).LE.epsilon)  GO TO 550
280         yo1      = yi
281         yi2      = yi * yi
282         Yf1      = a0 +  a1 * yi +     a2 * yi2  +     a3 * yi2 * yi
283         Yprimin  =       a1      + 2.* a2 *  yi  + 3.* a3 * yi2
284500   CONTINUE
285        WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter
286         STOP 2
287550   CONTINUE
288c
289       Yprimin   = a1  + 2.* a2 *  yi   + 3.* a3 * yi* yi
290       yprim(j)  = pi / ( jjm * Yprimin )
291       yvrai(j)  = yi
292
2931500    CONTINUE
294
295       DO j = 1, jlat -1
296        IF( yvrai(j+1). LT. yvrai(j) )  THEN
297         WRITE(6,*) ' PBS. avec  rlat(',j+1,') plus petit que rlat(',j,
298     ,  ')'
299         STOP 3
300        ENDIF
301       ENDDO
302
303       WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2'
304     , ,' et  pi/2 '
305c
306        IF( ik.EQ.1 )   THEN
307           ypn = pis2
308          DO j = jlat,1,-1
309           IF( yvrai(j).LE. ypn ) GO TO 1502
310          ENDDO
3111502     CONTINUE
312
313         jpn   = j
314         y00   = yvrai(jpn)
315         deply = pis2 -  y00
316        ENDIF
317
318         DO  j = 1, jjm +1 - jpn
319           ylatt (j)  = -pis2 - y00  + yvrai(jpn+j-1)
320           yprimm(j)  = yprim(jpn+j-1)
321         ENDDO
322
323         jjpn  = jpn
324         IF( jlat.EQ. jjm ) jjpn = jpn -1
325
326         DO j = 1,jjpn
327          ylatt (j + jjm+1 -jpn) = yvrai(j) + deply
328          yprimm(j + jjm+1 -jpn) = yprim(j)
329         ENDDO
330
331c      ***********   Fin de la reorganisation     *************
332c
333 1600   CONTINUE
334
335       DO j = 1, jlat
336          ylat(j) =  ylatt( jlat +1 -j )
337         yprim(j) = yprimm( jlat +1 -j )
338       ENDDO
339 
340        DO j = 1, jlat
341         yvrai(j) = ylat(j)*180./pi
342        ENDDO
343
344        IF( ik.EQ.1 )  THEN
345c         WRITE(6,18)
346c         WRITE(6,*)  ' YLAT  en U   apres ( en  deg. ) '
347c         WRITE(6,68) (yvrai(j),j=1,jlat)
348cc         WRITE(6,*) ' YPRIM '
349cc         WRITE(6,445) ( yprim(j),j=1,jlat)
350
351          DO j = 1, jlat
352            rrlatu(j) =  ylat( j )
353           yyprimu(j) = yprim( j )
354          ENDDO
355
356        ELSE IF ( ik.EQ. 2 )  THEN
357c         WRITE(6,18)
358c         WRITE(6,*) ' YLAT   en V  apres ( en  deg. ) '
359c         WRITE(6,68) (yvrai(j),j=1,jlat)
360cc         WRITE(6,*)' YPRIM '
361cc         WRITE(6,445) ( yprim(j),j=1,jlat)
362
363          DO j = 1, jlat
364            rrlatv(j) =  ylat( j )
365           yyprimv(j) = yprim( j )
366          ENDDO
367
368        ELSE IF ( ik.EQ. 3 )  THEN
369c         WRITE(6,18)
370c         WRITE(6,*)  ' YLAT  en U + 0.75  apres ( en  deg. ) '
371c         WRITE(6,68) (yvrai(j),j=1,jlat)
372cc         WRITE(6,*) ' YPRIM '
373cc         WRITE(6,445) ( yprim(j),j=1,jlat)
374
375          DO j = 1, jlat
376            rlatu2(j) =  ylat( j )
377           yprimu2(j) = yprim( j )
378          ENDDO
379
380        ELSE IF ( ik.EQ. 4 )  THEN
381c         WRITE(6,18)
382c         WRITE(6,*)  ' YLAT en U + 0.25  apres ( en  deg. ) '
383c         WRITE(6,68)(yvrai(j),j=1,jlat)
384cc         WRITE(6,*) ' YPRIM '
385cc         WRITE(6,68) ( yprim(j),j=1,jlat)
386
387          DO j = 1, jlat
388            rlatu1(j) =  ylat( j )
389           yprimu1(j) = yprim( j )
390          ENDDO
391
392        ENDIF
393
3945000   CONTINUE
395c
396        WRITE(6,18)
397c
398c  .....     fin de la boucle  do 5000 .....
399
400        DO j = 1, jjm
401         ylat(j) = rrlatu(j) - rrlatu(j+1)
402        ENDDO
403        champmin =  1.e12
404        champmax = -1.e12
405        DO j = 1, jjm
406         champmin = MIN( champmin, ylat(j) )
407         champmax = MAX( champmax, ylat(j) )
408        ENDDO
409         champmin = champmin * 180./pi
410         champmax = champmax * 180./pi
411
41224     FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3)
41318      FORMAT(/)
41468      FORMAT(1x,7f9.2)
415
416        RETURN
417        END
Note: See TracBrowser for help on using the repository browser.