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

Diff of /trunk/Sources/dyn3d/fyhyp.f

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

trunk/dyn3d/fyhyp.f revision 112 by guez, Thu Sep 18 13:36:51 2014 UTC trunk/Sources/dyn3d/fyhyp.f revision 149 by guez, Thu Jun 18 12:23:44 2015 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(jjm), yprimu2(jjm), rlatu1(jjm), 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        DOUBLE PRECISION champmin, champmax
30      INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax      INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax
31      REAL dzoom ! distance totale de la zone du zoom (en radians)      REAL dzoom ! distance totale de la zone du zoom (en radians)
32      DOUBLE PRECISION ylat(jjp1), yprim(jjp1)      DOUBLE PRECISION ylat(jjm + 1), yprim(jjm + 1)
33      DOUBLE PRECISION yuv      DOUBLE PRECISION yuv
34      DOUBLE PRECISION, save:: yt(0:nmax2)      DOUBLE PRECISION, save:: yt(0:nmax2)
35      DOUBLE PRECISION fhyp(0:nmax2), beta      DOUBLE PRECISION fhyp(0:nmax2), beta
# Line 48  contains Line 37  contains
37      DOUBLE PRECISION fxm(0:nmax2)      DOUBLE PRECISION fxm(0:nmax2)
38      DOUBLE PRECISION, save:: yf(0:nmax2)      DOUBLE PRECISION, save:: yf(0:nmax2)
39      DOUBLE PRECISION yypr(0:nmax2)      DOUBLE PRECISION yypr(0:nmax2)
40      DOUBLE PRECISION yvrai(jjp1), yprimm(jjp1), ylatt(jjp1)      DOUBLE PRECISION yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1)
41      DOUBLE PRECISION pi, pis2, epsilon, y0, pisjm      DOUBLE PRECISION pi, pis2, epsilon, pisjm
42      DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin      DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin
43      DOUBLE PRECISION yfi, yf1, ffdy      DOUBLE PRECISION yfi, yf1, ffdy
44      DOUBLE PRECISION ypn, deply, y00      DOUBLE PRECISION ypn, deply, y00
# Line 58  contains Line 47  contains
47      INTEGER i, j, it, ik, iter, jlat      INTEGER i, j, it, ik, iter, jlat
48      INTEGER jpn, jjpn      INTEGER jpn, jjpn
49      SAVE jpn      SAVE jpn
50      DOUBLE PRECISION a0, a1, a2, a3, yi2, heavyy0, heavyy0m      DOUBLE PRECISION yi2, heavyy0, heavyy0m
51      DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2)      DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2)
52      REAL y0min, y0max      REAL y0min, y0max
53    
     DOUBLE PRECISION heavyside  
   
54      !-------------------------------------------------------------------      !-------------------------------------------------------------------
55    
56        print *, "Call sequence information: fyhyp"
57    
58      pi = 2.*asin(1.)      pi = 2.*asin(1.)
59      pis2 = pi/2.      pis2 = pi/2.
60      pisjm = pi/real(jjm)      pisjm = pi/real(jjm)
61      epsilon = 1e-3      epsilon = 1e-3
62      y0 = yzoomdeg*pi/180.      dzoom = dzoomy*pi
63        print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'
64      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  
65    
66      DO i = 0, nmax2      DO i = 0, nmax2
67         yt(i) = -pis2 + real(i)*pi/nmax2         yt(i) = -pis2 + real(i)*pi/nmax2
68      END DO      END DO
69    
70      heavyy0m = heavyside(-y0)      heavyy0m = heavyside(-clat)
71      heavyy0 = heavyside(y0)      heavyy0 = heavyside(clat)
72      y0min = 2.*y0*heavyy0m - pis2      y0min = 2.*clat*heavyy0m - pis2
73      y0max = 2.*y0*heavyy0 + pis2      y0max = 2.*clat*heavyy0 + pis2
74    
75      fa = 999.999      fa = 999.999
76      fb = 999.999      fb = 999.999
77    
78      DO i = 0, nmax2      DO i = 0, nmax2
79         IF (yt(i)<y0) THEN         IF (yt(i)<clat) THEN
80            fa(i) = tau*(yt(i)-y0 + dzoom/2.)            fa(i) = tauy*(yt(i)-clat + dzoom/2.)
81            fb(i) = (yt(i)-2.*y0*heavyy0m + pis2)*(y0-yt(i))            fb(i) = (yt(i)-2.*clat*heavyy0m + pis2)*(clat-yt(i))
82         ELSE IF (yt(i)>y0) THEN         ELSE IF (yt(i)>clat) THEN
83            fa(i) = tau*(y0-yt(i) + dzoom/2.)            fa(i) = tauy*(clat-yt(i) + dzoom/2.)
84            fb(i) = (2.*y0*heavyy0-yt(i) + pis2)*(yt(i)-y0)            fb(i) = (2.*clat*heavyy0-yt(i) + pis2)*(yt(i)-clat)
85         END IF         END IF
86    
87         IF (200.*fb(i)<-fa(i)) THEN         IF (200.*fb(i)<-fa(i)) THEN
# Line 114  contains Line 92  contains
92            fhyp(i) = tanh(fa(i)/fb(i))            fhyp(i) = tanh(fa(i)/fb(i))
93         END IF         END IF
94    
95         IF (yt(i)==y0) fhyp(i) = 1.         IF (yt(i)==clat) fhyp(i) = 1.
96         IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.         IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1.
97      END DO      END DO
98    
# Line 124  contains Line 102  contains
102    
103      DO i = 1, nmax2      DO i = 1, nmax2
104         ymoy = 0.5*(yt(i-1) + yt(i))         ymoy = 0.5*(yt(i-1) + yt(i))
105         IF (ymoy<y0) THEN         IF (ymoy<clat) THEN
106            fa(i) = tau*(ymoy-y0 + dzoom/2.)            fa(i) = tauy*(ymoy-clat + dzoom/2.)
107            fb(i) = (ymoy-2.*y0*heavyy0m + pis2)*(y0-ymoy)            fb(i) = (ymoy-2.*clat*heavyy0m + pis2)*(clat-ymoy)
108         ELSE IF (ymoy>y0) THEN         ELSE IF (ymoy>clat) THEN
109            fa(i) = tau*(y0-ymoy + dzoom/2.)            fa(i) = tauy*(clat-ymoy + dzoom/2.)
110            fb(i) = (2.*y0*heavyy0-ymoy + pis2)*(ymoy-y0)            fb(i) = (2.*clat*heavyy0-ymoy + pis2)*(ymoy-clat)
111         END IF         END IF
112    
113         IF (200.*fb(i)<-fa(i)) THEN         IF (200.*fb(i)<-fa(i)) THEN
# Line 139  contains Line 117  contains
117         ELSE         ELSE
118            fxm(i) = tanh(fa(i)/fb(i))            fxm(i) = tanh(fa(i)/fb(i))
119         END IF         END IF
120         IF (ymoy==y0) fxm(i) = 1.         IF (ymoy==clat) fxm(i) = 1.
121         IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.         IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1.
122         ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))         ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1))
123      END DO      END DO
124    
125      beta = (grossism*ffdy-pi)/(ffdy-pi)      beta = (grossismy*ffdy-pi)/(ffdy-pi)
126    
127      IF (2. * beta - grossism <= 0.) THEN      IF (2. * beta - grossismy <= 0.) THEN
128         print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &         print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' &
129              // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &              // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' &
130              // 'dzoomy et relancer.'              // 'dzoomy et relancer.'
# Line 156  contains Line 134  contains
134      ! calcul de Ytprim      ! calcul de Ytprim
135    
136      DO i = 0, nmax2      DO i = 0, nmax2
137         ytprim(i) = beta + (grossism-beta)*fhyp(i)         ytprim(i) = beta + (grossismy-beta)*fhyp(i)
138      END DO      END DO
139    
140      ! Calcul de Yf      ! Calcul de Yf
141    
142      yf(0) = -pis2      yf(0) = -pis2
143      DO i = 1, nmax2      DO i = 1, nmax2
144         yypr(i) = beta + (grossism-beta)*fxm(i)         yypr(i) = beta + (grossismy-beta)*fxm(i)
145      END DO      END DO
146    
147      DO i = 1, nmax2      DO i = 1, nmax2
# Line 209  contains Line 187  contains
187            ! et Y'(yi)            ! et Y'(yi)
188    
189            CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &            CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), &
190                 yt(it), yt(it + 1), a0, a1, a2, a3)                 yt(it), yt(it + 1))
191    
192            yf1 = yf(it)            yf1 = yf(it)
193            yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi            yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi
# Line 245  contains Line 223  contains
223    
224         IF (ik==1) THEN         IF (ik==1) THEN
225            ypn = pis2            ypn = pis2
226            DO j = jlat, 1, -1            DO j = jjm + 1, 1, -1
227               IF (yvrai(j)<=ypn) exit               IF (yvrai(j)<=ypn) exit
228            END DO            END DO
229    
# Line 279  contains Line 257  contains
257         END DO         END DO
258    
259         IF (ik==1) THEN         IF (ik==1) THEN
260            DO j = 1, jlat            DO j = 1, jjm + 1
261               rrlatu(j) = ylat(j)               rlatu(j) = ylat(j)
              yyprimu(j) = yprim(j)  
262            END DO            END DO
263         ELSE IF (ik==2) THEN         ELSE IF (ik==2) THEN
264            DO j = 1, jlat            DO j = 1, jjm
265               rrlatv(j) = ylat(j)               rlatv(j) = ylat(j)
              yyprimv(j) = yprim(j)  
266            END DO            END DO
267         ELSE IF (ik==3) THEN         ELSE IF (ik==3) THEN
268            DO j = 1, jlat            DO j = 1, jjm
269               rlatu2(j) = ylat(j)               rlatu2(j) = ylat(j)
270               yprimu2(j) = yprim(j)               yprimu2(j) = yprim(j)
271            END DO            END DO
272         ELSE IF (ik==4) THEN         ELSE IF (ik==4) THEN
273            DO j = 1, jlat            DO j = 1, jjm
274               rlatu1(j) = ylat(j)               rlatu1(j) = ylat(j)
275               yprimu1(j) = yprim(j)               yprimu1(j) = yprim(j)
276            END DO            END DO
# Line 302  contains Line 278  contains
278      END DO loop_ik      END DO loop_ik
279    
280      DO j = 1, jjm      DO j = 1, jjm
281         ylat(j) = rrlatu(j) - rrlatu(j + 1)         ylat(j) = rlatu(j) - rlatu(j + 1)
282      END DO      END DO
283      champmin = 1e12      champmin = 1e12
284      champmax = -1e12      champmax = -1e12
# Line 313  contains Line 289  contains
289      champmin = champmin*180./pi      champmin = champmin*180./pi
290      champmax = champmax*180./pi      champmax = champmax*180./pi
291    
292        DO j = 1, jjm
293           IF (rlatu1(j) <= rlatu2(j)) THEN
294              print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j
295              STOP 13
296           ENDIF
297    
298           IF (rlatu2(j) <= rlatu(j+1)) THEN
299              print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j
300              STOP 14
301           ENDIF
302    
303           IF (rlatu(j) <= rlatu1(j)) THEN
304              print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j
305              STOP 15
306           ENDIF
307    
308           IF (rlatv(j) <= rlatu2(j)) THEN
309              print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j
310              STOP 16
311           ENDIF
312    
313           IF (rlatv(j) >= rlatu1(j)) THEN
314              print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j
315              STOP 17
316           ENDIF
317    
318           IF (rlatv(j) >= rlatu(j)) THEN
319              print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j
320              STOP 18
321           ENDIF
322        ENDDO
323    
324        print *, 'Latitudes'
325        print 3, champmin, champmax
326    
327    3   Format(1x, ' Au centre du zoom, la longueur de la maille est', &
328             ' d environ ', f0.2, ' degres ', /, &
329             ' alors que la maille en dehors de la zone du zoom est ', &
330             "d'environ ", f0.2, ' degres ')
331    
332    END SUBROUTINE fyhyp    END SUBROUTINE fyhyp
333    
334  end module fyhyp_m  end module fyhyp_m
Legend:
Removed from v.112  
changed lines
  Added in v.149
  ViewVC Help
Powered by ViewVC 1.1.21