source: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/fyhyp.F @ 5111

Last change on this file since 5111 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 KB
RevLine 
[524]1!
[1299]2! $Id: fyhyp.F 1299 2010-01-20 14:27:21Z abarral $
[524]3!
4c
5c
[650]6       SUBROUTINE fyhyp ( yzoomdeg, grossism, dzooma,tau  , 
[524]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#include "dimensions.h"
31#include "paramet.h"
32
33       INTEGER      nmax , nmax2
34       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
35c
36c
37c     .......  arguments  d'entree    .......
38c
[650]39       REAL yzoomdeg, grossism,dzooma,tau
[524]40c         ( rentres  par  run.def )
41
42c     .......  arguments  de sortie   .......
43c
44       REAL rrlatu(jjp1), yyprimu(jjp1),rrlatv(jjm), yyprimv(jjm),
45     , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)
46
47c
48c     .....     champs  locaux    .....
49c
50     
[650]51       REAL   dzoom
[1279]52       REAL(KIND=8) ylat(jjp1), yprim(jjp1)
53       REAL(KIND=8) yuv
54       REAL(KIND=8) yt(0:nmax2)
55       REAL(KIND=8) fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2)
[524]56       SAVE Ytprim, yt,Yf
[1279]57       REAL(KIND=8) Yf(0:nmax2),yypr(0:nmax2)
58       REAL(KIND=8) yvrai(jjp1), yprimm(jjp1),ylatt(jjp1)
59       REAL(KIND=8) pi,depi,pis2,epsilon,y0,pisjm
60       REAL(KIND=8) yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax
61       REAL(KIND=8) yfi,Yf1,ffdy
62       REAL(KIND=8) ypn,deply,y00
[524]63       SAVE y00, deply
64
65       INTEGER i,j,it,ik,iter,jlat
66       INTEGER jpn,jjpn
67       SAVE jpn
[1279]68       REAL(KIND=8) a0,a1,a2,a3,yi2,heavyy0,heavyy0m
69       REAL(KIND=8) fa(0:nmax2),fb(0:nmax2)
[524]70       REAL y0min,y0max
71
[1279]72       REAL(KIND=8)     heavyside
[524]73
74       pi       = 2. * ASIN(1.)
75       depi     = 2. * pi
76       pis2     = pi/2.
[1299]77       pisjm    = pi/ REAL(jjm)
[524]78       epsilon  = 1.e-3
79       y0       =  yzoomdeg * pi/180.
80
[650]81       IF( dzooma.LT.1.)  THEN
82         dzoom = dzooma * pi
83       ELSEIF( dzooma.LT. 12. ) THEN
[524]84         WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug
85     ,menter et relancer ! '
86         STOP 1
87       ELSE
[650]88         dzoom = dzooma * pi/180.
[524]89       ENDIF
90
91       WRITE(6,18)
92       WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)'
93       WRITE(6,24) y0,grossism,tau,dzoom
94
95       DO i = 0, nmax2
[1299]96        yt(i) = - pis2  + REAL(i)* pi /nmax2
[524]97       ENDDO
98
99       heavyy0m = heavyside( -y0 )
100       heavyy0  = heavyside(  y0 )
101       y0min    = 2.*y0*heavyy0m - pis2
102       y0max    = 2.*y0*heavyy0  + pis2
103
104       fa = 999.999
105       fb = 999.999
106       
107       DO i = 0, nmax2
108        IF( yt(i).LT.y0 )  THEN
109         fa (i) = tau*  (yt(i)-y0+dzoom/2. )
110         fb(i) =   (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) )
111        ELSEIF ( yt(i).GT.y0 )  THEN
112         fa(i) =   tau *(y0-yt(i)+dzoom/2. )
113         fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 )
114       ENDIF
115       
116       IF( 200.* fb(i) .LT. - fa(i) )   THEN
117         fhyp ( i) = - 1.
118       ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
119         fhyp ( i) =   1.
120       ELSE 
121         fhyp(i) =  TANH ( fa(i)/fb(i) )
122       ENDIF
123
124       IF( yt(i).EQ.y0 )  fhyp(i) = 1.
125       IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1.
126
127       ENDDO
128
129cc  ....  Calcul  de  beta  ....
130c
131       ffdy   = 0.
132
133       DO i = 1, nmax2
134        ymoy    = 0.5 * ( yt(i-1) + yt( i ) )
135        IF( ymoy.LT.y0 )  THEN
136         fa(i)= tau * ( ymoy-y0+dzoom/2.)
137         fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy )
138        ELSEIF ( ymoy.GT.y0 )  THEN
139         fa(i)= tau * ( y0-ymoy+dzoom/2. )
140         fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 )
141        ENDIF
142
143        IF( 200.* fb(i) .LT. - fa(i) )    THEN
144         fxm ( i) = - 1.
145        ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
146         fxm ( i) =   1.
147        ELSE
148         fxm(i) =  TANH ( fa(i)/fb(i) )
149        ENDIF
150         IF( ymoy.EQ.y0 )  fxm(i) = 1.
151         IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1.
152         ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) )
153
154        ENDDO
155
156        beta  = ( grossism * ffdy - pi ) / ( ffdy - pi )
157
158       IF( 2.*beta - grossism.LE. 0.)  THEN
159
160        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
161     ,tine fyhyp est mauvaise ! '
162        WRITE(6,*)'Modifier les valeurs de  grossismy ,tauy ou dzoomy',
163     , ' et relancer ! ***  '
164        CALL ABORT
165
166       ENDIF
167c
168c   .....  calcul  de  Ytprim   .....
169c
170       
171       DO i = 0, nmax2
172        Ytprim(i) = beta  + ( grossism - beta ) * fhyp(i)
173       ENDDO
174
175c   .....  Calcul  de  Yf     ........
176
177       Yf(0) = - pis2
178       DO i = 1, nmax2
179        yypr(i)    = beta + ( grossism - beta ) * fxm(i)
180       ENDDO
181
182       DO i=1,nmax2
183        Yf(i)   = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) )
184       ENDDO
185
186c    ****************************************************************
187c
188c   .....   yuv  = 0.   si calcul des latitudes  aux pts.  U  .....
189c   .....   yuv  = 0.5  si calcul des latitudes  aux pts.  V  .....
190c
191      WRITE(6,18)
192c
193      DO 5000  ik = 1,4
194
195       IF( ik.EQ.1 )  THEN
196         yuv  = 0.
197         jlat = jjm + 1
198       ELSE IF ( ik.EQ.2 )  THEN
199         yuv  = 0.5
200         jlat = jjm
201       ELSE IF ( ik.EQ.3 )  THEN
202         yuv  = 0.25
203         jlat = jjm
204       ELSE IF ( ik.EQ.4 )  THEN
205         yuv  = 0.75
206         jlat = jjm
207       ENDIF
208c
209       yo1   = 0.
210       DO 1500 j =  1,jlat
211        yo1   = 0.
[1299]212        ylon2 =  - pis2 + pisjm * ( REAL(j)  + yuv  -1.) 
[524]213        yfi    = ylon2
214c
215       DO 250 it =  nmax2,0,-1
216        IF( yfi.GE.Yf(it))  GO TO 350
217250    CONTINUE
218       it = 0
219350    CONTINUE
220
221       yi = yt(it)
222       IF(it.EQ.nmax2)  THEN
223        it       = nmax2 -1
224        Yf(it+1) = pis2
225       ENDIF
226c  .................................................................
227c  ....  Interpolation entre  yi(it) et yi(it+1)   pour avoir Y(yi) 
228c      .....           et   Y'(yi)                             .....
229c  .................................................................
230
231       CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1),   
232     ,                  yt(it),yt(it+1) ,   a0,a1,a2,a3   )     
233
234       Yf1     = Yf(it)
235       Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi
236
237       DO 500 iter = 1,300
238         yi = yi - ( Yf1 - yfi )/ Yprimin
239
240        IF( ABS(yi-yo1).LE.epsilon)  GO TO 550
241         yo1      = yi
242         yi2      = yi * yi
243         Yf1      = a0 +  a1 * yi +     a2 * yi2  +     a3 * yi2 * yi
244         Yprimin  =       a1      + 2.* a2 *  yi  + 3.* a3 * yi2
245500   CONTINUE
246        WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter
247         STOP 2
248550   CONTINUE
249c
250       Yprimin   = a1  + 2.* a2 *  yi   + 3.* a3 * yi* yi
251       yprim(j)  = pi / ( jjm * Yprimin )
252       yvrai(j)  = yi
253
2541500    CONTINUE
255
256       DO j = 1, jlat -1
257        IF( yvrai(j+1). LT. yvrai(j) )  THEN
258         WRITE(6,*) ' PBS. avec  rlat(',j+1,') plus petit que rlat(',j,
259     ,  ')'
260         STOP 3
261        ENDIF
262       ENDDO
263
264       WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2'
265     , ,' et  pi/2 '
266c
267        IF( ik.EQ.1 )   THEN
268           ypn = pis2
269          DO j = jlat,1,-1
270           IF( yvrai(j).LE. ypn ) GO TO 1502
271          ENDDO
2721502     CONTINUE
273
274         jpn   = j
275         y00   = yvrai(jpn)
276         deply = pis2 -  y00
277        ENDIF
278
279         DO  j = 1, jjm +1 - jpn
280           ylatt (j)  = -pis2 - y00  + yvrai(jpn+j-1)
281           yprimm(j)  = yprim(jpn+j-1)
282         ENDDO
283
284         jjpn  = jpn
285         IF( jlat.EQ. jjm ) jjpn = jpn -1
286
287         DO j = 1,jjpn
288          ylatt (j + jjm+1 -jpn) = yvrai(j) + deply
289          yprimm(j + jjm+1 -jpn) = yprim(j)
290         ENDDO
291
292c      ***********   Fin de la reorganisation     *************
293c
294 1600   CONTINUE
295
296       DO j = 1, jlat
297          ylat(j) =  ylatt( jlat +1 -j )
298         yprim(j) = yprimm( jlat +1 -j )
299       ENDDO
300 
301        DO j = 1, jlat
302         yvrai(j) = ylat(j)*180./pi
303        ENDDO
304
305        IF( ik.EQ.1 )  THEN
306c         WRITE(6,18)
307c         WRITE(6,*)  ' YLAT  en U   apres ( en  deg. ) '
308c         WRITE(6,68) (yvrai(j),j=1,jlat)
309cc         WRITE(6,*) ' YPRIM '
310cc         WRITE(6,445) ( yprim(j),j=1,jlat)
311
312          DO j = 1, jlat
313            rrlatu(j) =  ylat( j )
314           yyprimu(j) = yprim( j )
315          ENDDO
316
317        ELSE IF ( ik.EQ. 2 )  THEN
318c         WRITE(6,18)
319c         WRITE(6,*) ' YLAT   en V  apres ( en  deg. ) '
320c         WRITE(6,68) (yvrai(j),j=1,jlat)
321cc         WRITE(6,*)' YPRIM '
322cc         WRITE(6,445) ( yprim(j),j=1,jlat)
323
324          DO j = 1, jlat
325            rrlatv(j) =  ylat( j )
326           yyprimv(j) = yprim( j )
327          ENDDO
328
329        ELSE IF ( ik.EQ. 3 )  THEN
330c         WRITE(6,18)
331c         WRITE(6,*)  ' YLAT  en U + 0.75  apres ( en  deg. ) '
332c         WRITE(6,68) (yvrai(j),j=1,jlat)
333cc         WRITE(6,*) ' YPRIM '
334cc         WRITE(6,445) ( yprim(j),j=1,jlat)
335
336          DO j = 1, jlat
337            rlatu2(j) =  ylat( j )
338           yprimu2(j) = yprim( j )
339          ENDDO
340
341        ELSE IF ( ik.EQ. 4 )  THEN
342c         WRITE(6,18)
343c         WRITE(6,*)  ' YLAT en U + 0.25  apres ( en  deg. ) '
344c         WRITE(6,68)(yvrai(j),j=1,jlat)
345cc         WRITE(6,*) ' YPRIM '
346cc         WRITE(6,68) ( yprim(j),j=1,jlat)
347
348          DO j = 1, jlat
349            rlatu1(j) =  ylat( j )
350           yprimu1(j) = yprim( j )
351          ENDDO
352
353        ENDIF
354
3555000   CONTINUE
356c
357        WRITE(6,18)
358c
359c  .....     fin de la boucle  do 5000 .....
360
361        DO j = 1, jjm
362         ylat(j) = rrlatu(j) - rrlatu(j+1)
363        ENDDO
364        champmin =  1.e12
365        champmax = -1.e12
366        DO j = 1, jjm
367         champmin = MIN( champmin, ylat(j) )
368         champmax = MAX( champmax, ylat(j) )
369        ENDDO
370         champmin = champmin * 180./pi
371         champmax = champmax * 180./pi
372
37324     FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3)
37418      FORMAT(/)
37568      FORMAT(1x,7f9.2)
376
377        RETURN
378        END
Note: See TracBrowser for help on using the repository browser.