/[lmdze]/trunk/dyn3d/fyhyp.f
ViewVC logotype

Diff of /trunk/dyn3d/fyhyp.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 112 by guez, Thu Sep 18 13:36:51 2014 UTC revision 254 by guez, Mon Feb 5 10:39:38 2018 UTC
# Line 4  module fyhyp_m Line 4  module fyhyp_m
4    
5  contains  contains
6    
7    SUBROUTINE fyhyp(yzoomdeg, grossism, dzooma, tau, rrlatu, yyprimu, rrlatv, &    SUBROUTINE fyhyp(rlatu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
        yyprimv, rlatu2, yprimu2, rlatu1, yprimu1, champmin, champmax)  
8    
9      ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32      ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32
10    
11      ! Author: P. Le Van, from analysis by R. Sadourny      ! Author: P. Le Van, from analysis by R. Sadourny
12    
13      ! Calcule les latitudes et dérivées dans la grille du GCM pour une      ! Calcule les latitudes et dérivées dans la grille du GCM pour une
14      ! fonction f(y) à tangente hyperbolique.      ! fonction f(y) à dérivée tangente hyperbolique.
15    
16      ! Nota bene : il vaut mieux avoir grossism * dzoom < pi / 2 (rad),      ! Il vaut mieux avoir : grossismy * dzoom < pi / 2
     ! en latitude.  
17    
18        use coefpoly_m, only: coefpoly, a0, a1, a2, a3
19      USE dimens_m, only: jjm      USE dimens_m, only: jjm
20      USE paramet_m, only: JJP1      use dynetat0_m, only: clat, grossismy, dzoomy, tauy
21        use heavyside_m, only: heavyside
22    
23      REAL, intent(in):: yzoomdeg      REAL, intent(out):: rlatu(:) ! (jjm + 1)
24        REAL, intent(out):: rlatv(:) ! (jjm)
25      REAL, intent(in):: grossism      real, intent(out):: rlatu2(:), yprimu2(:), rlatu1(:), yprimu1(:) ! (jjm)
     ! grossissement (= 2 si 2 fois, = 3 si 3 fois, etc.)  
   
     REAL, intent(in):: dzooma  
   
     REAL, intent(in):: tau  
     ! raideur de la transition de l'intérieur à l'extérieur du zoom  
   
     ! arguments de sortie  
   
     REAL rrlatu(jjp1), yyprimu(jjp1), rrlatv(jjm), yyprimv(jjm)  
     real rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm)  
     DOUBLE PRECISION champmin, champmax  
26    
27      ! Local:      ! Local:
28    
29      INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax      INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax
30      REAL dzoom ! distance totale de la zone du zoom (en radians)      REAL dzoom ! distance totale de la zone du zoom (en radians)
31      DOUBLE PRECISION ylat(jjp1), yprim(jjp1)      DOUBLE PRECISION ylat(jjm + 1), yprim(jjm + 1)
32      DOUBLE PRECISION yuv      DOUBLE PRECISION yuv
33      DOUBLE PRECISION, save:: yt(0:nmax2)      DOUBLE PRECISION, save:: yt(0:nmax2)
34      DOUBLE PRECISION fhyp(0:nmax2), beta      DOUBLE PRECISION fhyp(0:nmax2), beta
# Line 48  contains Line 36  contains
36      DOUBLE PRECISION fxm(0:nmax2)      DOUBLE PRECISION fxm(0:nmax2)
37      DOUBLE PRECISION, save:: yf(0:nmax2)      DOUBLE PRECISION, save:: yf(0:nmax2)
38      DOUBLE PRECISION yypr(0:nmax2)      DOUBLE PRECISION yypr(0:nmax2)
39      DOUBLE PRECISION yvrai(jjp1), yprimm(jjp1), ylatt(jjp1)      DOUBLE PRECISION yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1)
40      DOUBLE PRECISION pi, pis2, epsilon, y0, pisjm      DOUBLE PRECISION pi, pis2, epsilon, pisjm
41      DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin      DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin
42      DOUBLE PRECISION yfi, yf1, ffdy      DOUBLE PRECISION yfi, yf1, ffdy
43      DOUBLE PRECISION ypn, deply, y00      DOUBLE PRECISION ypn
44      SAVE y00, deply      DOUBLE PRECISION, save::deply, y00
45    
46      INTEGER i, j, it, ik, iter, jlat      INTEGER i, j, it, ik, iter, jlat, jjpn
47      INTEGER jpn, jjpn      INTEGER, save:: jpn
48      SAVE jpn      DOUBLE PRECISION yi2, heavyy0, heavyy0m
     DOUBLE PRECISION a0, a1, a2, a3, yi2, heavyy0, heavyy0m  
49      DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2)      DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2)
50      REAL y0min, y0max      REAL y0min, y0max
51    
     DOUBLE PRECISION heavyside  
   
52      !-------------------------------------------------------------------      !-------------------------------------------------------------------
53    
54        print *, "Call sequence information: fyhyp"
55    
56      pi = 2.*asin(1.)      pi = 2.*asin(1.)
57      pis2 = pi/2.      pis2 = pi/2.
58      pisjm = pi/real(jjm)      pisjm = pi/real(jjm)
59      epsilon = 1e-3      epsilon = 1e-3
60      y0 = yzoomdeg*pi/180.      dzoom = dzoomy*pi
61        print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'
62      IF (dzooma<1.) THEN      print *, clat, grossismy, tauy, dzoom
        dzoom = dzooma*pi  
     ELSE IF (dzooma<12.) THEN  
        print *, "Le paramètre dzoomy pour fyhyp est trop petit. L'augmenter " &  
             // "et relancer."  
        STOP 1  
     ELSE  
        dzoom = dzooma * pi/180.  
     END IF  
   
     print *, 'yzoom(rad), grossism, tau, dzoom (rad):'  
     print *, y0, grossism, tau, dzoom  
63    
64      DO i = 0, nmax2      DO i = 0, nmax2
65         yt(i) = -pis2 + real(i)*pi/nmax2         yt(i) = -pis2 + real(i)*pi/nmax2
66      END DO      END DO
67    
68      heavyy0m = heavyside(-y0)      heavyy0m = heavyside(-clat)
69      heavyy0 = heavyside(y0)      heavyy0 = heavyside(clat)
70      y0min = 2.*y0*heavyy0m - pis2      y0min = 2.*clat*heavyy0m - pis2
71      y0max = 2.*y0*heavyy0 + pis2      y0max = 2.*clat*heavyy0 + pis2
72    
73      fa = 999.999      fa = 999.999
74      fb = 999.999      fb = 999.999
75    
76      DO i = 0, nmax2      DO i = 0, nmax2
77         IF (yt(i)<y0) THEN         IF (yt(i)<clat) THEN
78            fa(i) = tau*(yt(i)-y0 + dzoom/2.)            fa(i) = tauy*(yt(i)-clat + dzoom/2.)
79            fb(i) = (yt(i)-2.*y0*heavyy0m + pis2)*(y0-yt(i))            fb(i) = (yt(i)-2.*clat*heavyy0m + pis2)*(clat-yt(i))
80         ELSE IF (yt(i)>y0) THEN         ELSE IF (yt(i)>clat) THEN
81            fa(i) = tau*(y0-yt(i) + dzoom/2.)            fa(i) = tauy*(clat-yt(i) + dzoom/2.)
82            fb(i) = (2.*y0*heavyy0-yt(i) + pis2)*(yt(i)-y0)            fb(i) = (2.*clat*heavyy0-yt(i) + pis2)*(yt(i)-clat)
83         END IF         END IF
84    
85         IF (200.*fb(i)<-fa(i)) THEN         IF (200.*fb(i)<-fa(i)) THEN
# Line 114  contains Line 90  contains
90            fhyp(i) = tanh(fa(i)/fb(i))            fhyp(i) = tanh(fa(i)/fb(i))
91         END IF         END IF
92    
93         IF (yt(i)==y0) fhyp(i) = 1.         IF (yt(i)==clat) fhyp(i) = 1.
94         IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.         IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.
95      END DO      END DO
96    
# Line 124  contains Line 100  contains
100    
101      DO i = 1, nmax2      DO i = 1, nmax2
102         ymoy = 0.5*(yt(i-1) + yt(i))         ymoy = 0.5*(yt(i-1) + yt(i))
103         IF (ymoy<y0) THEN         IF (ymoy<clat) THEN
104            fa(i) = tau*(ymoy-y0 + dzoom/2.)            fa(i) = tauy*(ymoy-clat + dzoom/2.)
105            fb(i) = (ymoy-2.*y0*heavyy0m + pis2)*(y0-ymoy)            fb(i) = (ymoy-2.*clat*heavyy0m + pis2)*(clat-ymoy)
106         ELSE IF (ymoy>y0) THEN         ELSE IF (ymoy>clat) THEN
107            fa(i) = tau*(y0-ymoy + dzoom/2.)            fa(i) = tauy*(clat-ymoy + dzoom/2.)
108            fb(i) = (2.*y0*heavyy0-ymoy + pis2)*(ymoy-y0)            fb(i) = (2.*clat*heavyy0-ymoy + pis2)*(ymoy-clat)
109         END IF         END IF
110    
111         IF (200.*fb(i)<-fa(i)) THEN         IF (200.*fb(i)<-fa(i)) THEN
# Line 139  contains Line 115  contains
115         ELSE         ELSE
116            fxm(i) = tanh(fa(i)/fb(i))            fxm(i) = tanh(fa(i)/fb(i))
117         END IF         END IF
118         IF (ymoy==y0) fxm(i) = 1.         IF (ymoy==clat) fxm(i) = 1.
119         IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.         IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.
120         ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))         ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))
121      END DO      END DO
122    
123      beta = (grossism*ffdy-pi)/(ffdy-pi)      beta = (grossismy*ffdy-pi)/(ffdy-pi)
124    
125      IF (2. * beta - grossism <= 0.) THEN      IF (2. * beta - grossismy <= 0.) THEN
126         print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &         print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &
127              // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &              // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &
128              // 'dzoomy et relancer.'              // 'dzoomy et relancer.'
# Line 156  contains Line 132  contains
132      ! calcul de Ytprim      ! calcul de Ytprim
133    
134      DO i = 0, nmax2      DO i = 0, nmax2
135         ytprim(i) = beta + (grossism-beta)*fhyp(i)         ytprim(i) = beta + (grossismy-beta)*fhyp(i)
136      END DO      END DO
137    
138      ! Calcul de Yf      ! Calcul de Yf
139    
140      yf(0) = -pis2      yf(0) = -pis2
141      DO i = 1, nmax2      DO i = 1, nmax2
142         yypr(i) = beta + (grossism-beta)*fxm(i)         yypr(i) = beta + (grossismy-beta)*fxm(i)
143      END DO      END DO
144    
145      DO i = 1, nmax2      DO i = 1, nmax2
# Line 209  contains Line 185  contains
185            ! et Y'(yi)            ! et Y'(yi)
186    
187            CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &            CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &
188                 yt(it), yt(it + 1), a0, a1, a2, a3)                 yt(it), yt(it + 1))
189    
190            yf1 = yf(it)            yf1 = yf(it)
191            yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi            yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi
# Line 245  contains Line 221  contains
221    
222         IF (ik==1) THEN         IF (ik==1) THEN
223            ypn = pis2            ypn = pis2
224            DO j = jlat, 1, -1            DO j = jjm + 1, 1, -1
225               IF (yvrai(j)<=ypn) exit               IF (yvrai(j)<=ypn) exit
226            END DO            END DO
227    
# Line 279  contains Line 255  contains
255         END DO         END DO
256    
257         IF (ik==1) THEN         IF (ik==1) THEN
258            DO j = 1, jlat            DO j = 1, jjm + 1
259               rrlatu(j) = ylat(j)               rlatu(j) = ylat(j)
              yyprimu(j) = yprim(j)  
260            END DO            END DO
261         ELSE IF (ik==2) THEN         ELSE IF (ik==2) THEN
262            DO j = 1, jlat            DO j = 1, jjm
263               rrlatv(j) = ylat(j)               rlatv(j) = ylat(j)
              yyprimv(j) = yprim(j)  
264            END DO            END DO
265         ELSE IF (ik==3) THEN         ELSE IF (ik==3) THEN
266            DO j = 1, jlat            DO j = 1, jjm
267               rlatu2(j) = ylat(j)               rlatu2(j) = ylat(j)
268               yprimu2(j) = yprim(j)               yprimu2(j) = yprim(j)
269            END DO            END DO
270         ELSE IF (ik==4) THEN         ELSE IF (ik==4) THEN
271            DO j = 1, jlat            DO j = 1, jjm
272               rlatu1(j) = ylat(j)               rlatu1(j) = ylat(j)
273               yprimu1(j) = yprim(j)               yprimu1(j) = yprim(j)
274            END DO            END DO
# Line 302  contains Line 276  contains
276      END DO loop_ik      END DO loop_ik
277    
278      DO j = 1, jjm      DO j = 1, jjm
279         ylat(j) = rrlatu(j) - rrlatu(j + 1)         ylat(j) = rlatu(j) - rlatu(j + 1)
280      END DO      END DO
281      champmin = 1e12  
     champmax = -1e12  
282      DO j = 1, jjm      DO j = 1, jjm
283         champmin = min(champmin, ylat(j))         IF (rlatu1(j) <= rlatu2(j)) THEN
284         champmax = max(champmax, ylat(j))            print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j
285      END DO            STOP 13
286      champmin = champmin*180./pi         ENDIF
287      champmax = champmax*180./pi  
288           IF (rlatu2(j) <= rlatu(j+1)) THEN
289              print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j
290              STOP 14
291           ENDIF
292    
293           IF (rlatu(j) <= rlatu1(j)) THEN
294              print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j
295              STOP 15
296           ENDIF
297    
298           IF (rlatv(j) <= rlatu2(j)) THEN
299              print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j
300              STOP 16
301           ENDIF
302    
303           IF (rlatv(j) >= rlatu1(j)) THEN
304              print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j
305              STOP 17
306           ENDIF
307    
308           IF (rlatv(j) >= rlatu(j)) THEN
309              print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j
310              STOP 18
311           ENDIF
312        ENDDO
313    
314        print *, 'Latitudes'
315        print 3, minval(ylat(:jjm)) *180d0/pi, maxval(ylat(:jjm))*180d0/pi
316    
317    3   Format(1x, ' Au centre du zoom, la longueur de la maille est', &
318             ' d environ ', f0.2, ' degres ', /, &
319             ' alors que la maille en dehors de la zone du zoom est ', &
320             "d'environ ", f0.2, ' degres ')
321    
322    END SUBROUTINE fyhyp    END SUBROUTINE fyhyp
323    
Legend:
Removed from v.112  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21