/[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 121 by guez, Wed Jan 28 16:10:02 2015 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(rlatu, yyprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)    SUBROUTINE fyhyp(rlatu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
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    
# Line 15  contains Line 15  contains
15    
16      ! Il vaut mieux avoir : grossismy * dzoom < pi / 2      ! Il vaut mieux avoir : grossismy * dzoom < pi / 2
17    
18      use coefpoly_m, only: coefpoly      use coefpoly_m, only: coefpoly, a0, a1, a2, a3
19      USE dimens_m, only: jjm      USE dimens_m, only: jjm
20      use serre, only: clat, grossismy, dzoomy, tauy      use dynetat0_m, only: clat, grossismy, dzoomy, tauy
21        use heavyside_m, only: heavyside
22    
23      REAL, intent(out):: rlatu(jjm + 1), yyprimu(jjm + 1)      REAL, intent(out):: rlatu(:) ! (jjm + 1)
24      REAL, intent(out):: rlatv(jjm)      REAL, intent(out):: rlatv(:) ! (jjm)
25      real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm)      real, intent(out):: rlatu2(:), yprimu2(:), rlatu1(:), yprimu1(:) ! (jjm)
26    
27      ! Local:      ! Local:
28    
     DOUBLE PRECISION champmin, champmax  
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(jjm + 1), yprim(jjm + 1)      DOUBLE PRECISION ylat(jjm + 1), yprim(jjm + 1)
# Line 37  contains Line 37  contains
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(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1)      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"      print *, "Call sequence information: fyhyp"
# Line 60  contains Line 57  contains
57      pis2 = pi/2.      pis2 = pi/2.
58      pisjm = pi/real(jjm)      pisjm = pi/real(jjm)
59      epsilon = 1e-3      epsilon = 1e-3
     y0 = clat*pi/180.  
60      dzoom = dzoomy*pi      dzoom = dzoomy*pi
61      print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'      print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):'
62      print *, y0, grossismy, tauy, dzoom      print *, clat, grossismy, tauy, 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) = tauy*(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) = tauy*(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 94  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 104  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) = tauy*(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) = tauy*(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 119  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
# Line 189  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 261  contains Line 257  contains
257         IF (ik==1) THEN         IF (ik==1) THEN
258            DO j = 1, jjm + 1            DO j = 1, jjm + 1
259               rlatu(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, jjm            DO j = 1, jjm
# Line 283  contains Line 278  contains
278      DO j = 1, jjm      DO j = 1, jjm
279         ylat(j) = rlatu(j) - rlatu(j + 1)         ylat(j) = rlatu(j) - rlatu(j + 1)
280      END DO      END DO
     champmin = 1e12  
     champmax = -1e12  
     DO j = 1, jjm  
        champmin = min(champmin, ylat(j))  
        champmax = max(champmax, ylat(j))  
     END DO  
     champmin = champmin*180./pi  
     champmax = champmax*180./pi  
281    
282      DO j = 1, jjm      DO j = 1, jjm
283         IF (rlatu1(j) <= rlatu2(j)) THEN         IF (rlatu1(j) <= rlatu2(j)) THEN
# Line 325  contains Line 312  contains
312      ENDDO      ENDDO
313    
314      print *, 'Latitudes'      print *, 'Latitudes'
315      print 3, champmin, champmax      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', &  3   Format(1x, ' Au centre du zoom, la longueur de la maille est', &
318           ' d environ ', f0.2, ' degres ', /, &           ' d environ ', f0.2, ' degres ', /, &
Legend:
Removed from v.121  
changed lines
  Added in v.254
  ViewVC Help
Powered by ViewVC 1.1.21