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

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

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

trunk/dyn3d/comgeom.f revision 96 by guez, Fri Apr 4 11:30:34 2014 UTC trunk/Sources/dyn3d/comgeom.f revision 178 by guez, Fri Mar 11 18:47:26 2016 UTC
# Line 64  module comgeom Line 64  module comgeom
64    real fext((iim + 1) * jjm), constang((iim + 1) * (jjm + 1))    real fext((iim + 1) * jjm), constang((iim + 1) * (jjm + 1))
65    equivalence (fext, fext_2d), (constang, constang_2d)    equivalence (fext, fext_2d), (constang, constang_2d)
66    
   real rlatu(jjm + 1)  
   ! (latitudes of points of the "scalar" and "u" grid, in rad)  
   
   real rlatv(jjm)  
   ! (latitudes of points of the "v" grid, in rad, in decreasing order)  
   
   real rlonu(iim + 1) ! longitudes of points of the "u" grid, in rad  
   
   real rlonv(iim + 1)  
   ! (longitudes of points of the "scalar" and "v" grid, in rad)  
   
67    real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension    real cuvsurcv_2d(iim + 1, jjm), cvsurcuv_2d(iim + 1, jjm) ! no dimension
68    real cuvsurcv((iim + 1) * jjm), cvsurcuv((iim + 1) * jjm) ! no dimension    real cuvsurcv((iim + 1) * jjm), cvsurcuv((iim + 1) * jjm) ! no dimension
69    equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)    equivalence (cuvsurcv, cuvsurcv_2d), (cvsurcuv, cvsurcuv_2d)
# Line 111  module comgeom Line 100  module comgeom
100    real unsairz_gam((iim + 1) * jjm)    real unsairz_gam((iim + 1) * jjm)
101    equivalence (unsairz_gam, unsairz_gam_2d)    equivalence (unsairz_gam, unsairz_gam_2d)
102    
   real xprimu(iim + 1), xprimv(iim + 1)  
   
103    save    save
104    
105  contains  contains
# Line 124  contains Line 111  contains
111      ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes      ! Calcul des élongations cuij1, ..., cuij4, cvij1, ..., cvij4 aux mêmes
112      ! endroits que les aires aireij1_2d, ..., aireij4_2d.      ! endroits que les aires aireij1_2d, ..., aireij4_2d.
113    
114      ! Choix entre une fonction "f(y)" à dérivée sinusoïdale ou à      ! Calcul des coefficients cu_2d, cv_2d, 1. / cu_2d**2, 1. /
115      ! dérivée tangente hyperbolique. Calcul des coefficients cu_2d,      ! cv_2d**2. Les coefficients cu_2d et cv_2d permettent de passer
116      ! cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les coefficients cu_2d et cv_2d      ! des vitesses naturelles aux vitesses covariantes et
117      ! permettent de passer des vitesses naturelles aux vitesses      ! contravariantes, ou vice-versa.
     ! covariantes et contravariantes, ou vice-versa.  
118    
119      ! On a :      ! On a :
120      ! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d      ! u(covariant) = cu_2d * u(naturel), u(contravariant) = u(naturel) / cu_2d
# Line 138  contains Line 124  contains
124      ! u(covariant) = cu_2d * cu_2d * u(contravariant)      ! u(covariant) = cu_2d * cu_2d * u(contravariant)
125      ! v(covariant) = cv_2d * cv_2d * v(contravariant)      ! v(covariant) = cv_2d * cv_2d * v(contravariant)
126    
     ! On a l'application (x(X), y(Y)) avec - im / 2 + 1 <= X <= im / 2  
     ! et - jm / 2 <= Y <= jm / 2  
   
127      ! x est la longitude du point en radians.      ! x est la longitude du point en radians.
128      ! y est la latitude du point en radians.      ! y est la latitude du point en radians.
129      !      !
# Line 152  contains Line 135  contains
135      ! dépendant de j uniquement, sera ici indicé aussi en i pour un      ! dépendant de j uniquement, sera ici indicé aussi en i pour un
136      ! adressage plus facile en ij.      ! adressage plus facile en ij.
137    
138      ! xprimu et xprimv sont respectivement les valeurs de dx / dX aux      ! cv_2d est aux points v. cu_2d est aux points u. Cf. "inigeom.txt".
     ! points u et v. yprimu et yprimv sont respectivement les valeurs  
     ! de dy / dY aux points u et v. rlatu et rlatv sont respectivement  
     ! les valeurs de la latitude aux points u et v. cvu et cv_2d sont  
     ! respectivement les valeurs de cv_2d aux points u et v.  
   
     ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d  
     ! aux points u, v, scalaires, et z. Cf. "inigeom.txt".  
139    
140      USE comconst, ONLY : g, omeg, rad      USE comconst, ONLY : g, omeg, rad
141      USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh      USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
142      use conf_gcm_m, ONLY : fxyhypb, ysinus      use dynetat0_m, only: xprimp025, xprimm025, rlatu1, rlatu2, rlatu, rlatv, &
143      use fxy_m, only: fxy           yprimu1, yprimu2
     use fxyhyper_m, only: fxyhyper  
     use jumble, only: new_unit  
     use nr_util, only: pi  
144      USE paramet_m, ONLY : iip1, jjp1      USE paramet_m, ONLY : iip1, jjp1
     USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, &  
          grossismy, pxo, pyo, taux, tauy, transx, transy  
     ! Modifiés pxo, pyo, transx, transy  
145    
146      ! Local:      ! Local:
147      INTEGER i, j, itmax, itmay, iter, unit      INTEGER i, j
     REAL cvu(iip1, jjp1), cuv(iip1, jjm)  
148      REAL ai14, ai23, airez, un4rad2      REAL ai14, ai23, airez, un4rad2
     REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm  
149      REAL coslatm, coslatp, radclatm, radclatp      REAL coslatm, coslatp, radclatm, radclatp
150      REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m      REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
151      REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m      REAL, dimension(iip1, jjp1):: cvij1, cvij2, cvij3, cvij4 ! in m
     REAL rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm)  
     real yprimv(jjm), yprimu(jjp1)  
152      REAL gamdi_gdiv, gamdi_grot, gamdi_h      REAL gamdi_gdiv, gamdi_grot, gamdi_h
     REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)  
153      real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &      real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
154           aireij4_2d ! in m2           aireij4_2d ! in m2
     real airuscv2_2d(iim + 1, jjm)  
     real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm)  
     real aivscu2gam_2d(iim + 1, jjm)  
155    
156      !------------------------------------------------------------------      !------------------------------------------------------------------
157    
158      PRINT *, 'Call sequence information: inigeom'      PRINT *, 'Call sequence information: inigeom'
159    
160      IF (nitergdiv/=2) THEN      IF (nitergdiv /= 2) THEN
161         gamdi_gdiv = coefdis / (real(nitergdiv)-2.)         gamdi_gdiv = coefdis / (nitergdiv - 2)
162      ELSE      ELSE
163         gamdi_gdiv = 0.         gamdi_gdiv = 0.
164      END IF      END IF
165      IF (nitergrot/=2) THEN  
166         gamdi_grot = coefdis / (real(nitergrot)-2.)      IF (nitergrot /= 2) THEN
167           gamdi_grot = coefdis / (nitergrot - 2)
168      ELSE      ELSE
169         gamdi_grot = 0.         gamdi_grot = 0.
170      END IF      END IF
171      IF (niterh/=2) THEN  
172         gamdi_h = coefdis / (real(niterh)-2.)      IF (niterh /= 2) THEN
173           gamdi_h = coefdis / (niterh - 2)
174      ELSE      ELSE
175         gamdi_h = 0.         gamdi_h = 0.
176      END IF      END IF
# Line 215  contains Line 179  contains
179      print *, "gamdi_grot = ", gamdi_grot      print *, "gamdi_grot = ", gamdi_grot
180      print *, "gamdi_h = ", gamdi_h      print *, "gamdi_h = ", gamdi_h
181    
     IF (.NOT. fxyhypb) THEN  
        IF (ysinus) THEN  
           print *, ' Inigeom, Y = Sinus (Latitude) '  
           ! utilisation de f(x, y) avec y = sinus de la latitude  
           CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, &  
                rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, &  
                xprimm025, rlonp025, xprimp025)  
        ELSE  
           print *, 'Inigeom, Y = Latitude, der. sinusoid .'  
           ! utilisation de f(x, y) a tangente sinusoidale, y etant la latit  
   
           pxo = clon * pi / 180.  
           pyo = 2. * clat * pi / 180.  
   
           ! determination de transx (pour le zoom) par Newton-Raphson  
   
           itmax = 10  
           eps = .1E-7  
   
           xo1 = 0.  
           DO iter = 1, itmax  
              x1 = xo1  
              f = x1 + alphax * sin(x1-pxo)  
              df = 1. + alphax * cos(x1-pxo)  
              x1 = x1 - f / df  
              xdm = abs(x1-xo1)  
              IF (xdm<=eps) EXIT  
              xo1 = x1  
           END DO  
   
           transx = xo1  
   
           itmay = 10  
           eps = .1E-7  
   
           yo1 = 0.  
           DO iter = 1, itmay  
              y1 = yo1  
              f = y1 + alphay * sin(y1-pyo)  
              df = 1. + alphay * cos(y1-pyo)  
              y1 = y1 - f / df  
              ydm = abs(y1-yo1)  
              IF (ydm<=eps) EXIT  
              yo1 = y1  
           END DO  
   
           transy = yo1  
   
           CALL fxy(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &  
                yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &  
                rlonp025, xprimp025)  
        END IF  
     ELSE  
        ! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique  
        print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique'  
        CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, &  
             taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, &  
             yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, &  
             rlonp025, xprimp025)  
     END IF  
   
     rlatu(1) = pi / 2.  
     rlatu(jjp1) = -rlatu(1)  
   
     ! Calcul aux pôles  
   
     yprimu(1) = 0.  
     yprimu(jjp1) = 0.  
   
182      un4rad2 = 0.25 * rad * rad      un4rad2 = 0.25 * rad * rad
183    
184      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
# Line 426  contains Line 321  contains
321         unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)         unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
322      END DO      END DO
323    
324      ! Calcul des élongations cu_2d, cv_2d, cvu      ! Calcul des élongations cu_2d, cv_2d
325    
326      DO j = 1, jjm      DO j = 1, jjm
327         DO i = 1, iim         DO i = 1, iim
328            cv_2d(i, j) = 0.5 * &            cv_2d(i, j) = 0.5 * &
329                 (cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))                 (cvij2(i, j) + cvij3(i, j) + cvij1(i, j + 1) + cvij4(i, j + 1))
           cvu(i, j) = 0.5 * (cvij1(i, j) + cvij4(i, j) + cvij2(i, j) &  
                + cvij3(i, j))  
           cuv(i, j) = 0.5 * (cuij2(i, j) + cuij3(i, j) + cuij1(i, j + 1) &  
                + cuij4(i, j + 1))  
330            unscv2_2d(i, j) = 1. / cv_2d(i, j)**2            unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
331         END DO         END DO
332         DO i = 1, iim         DO i = 1, iim
# Line 446  contains Line 337  contains
337            cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)            cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
338         END DO         END DO
339         cv_2d(iip1, j) = cv_2d(1, j)         cv_2d(iip1, j) = cv_2d(1, j)
        cvu(iip1, j) = cvu(1, j)  
340         unscv2_2d(iip1, j) = unscv2_2d(1, j)         unscv2_2d(iip1, j) = unscv2_2d(1, j)
        cuv(iip1, j) = cuv(1, j)  
341         cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)         cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
342         cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)         cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
343         cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)         cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
# Line 480  contains Line 369  contains
369    
370      cu_2d(:, 1) = 0.      cu_2d(:, 1) = 0.
371      unscu2_2d(:, 1) = 0.      unscu2_2d(:, 1) = 0.
     cvu(:, 1) = 0.  
372    
373      cu_2d(:, jjp1) = 0.      cu_2d(:, jjp1) = 0.
374      unscu2_2d(:, jjp1) = 0.      unscu2_2d(:, jjp1) = 0.
     cvu(:, jjp1) = 0.  
   
     DO j = 1, jjm  
        DO i = 1, iim  
           airvscu2_2d(i, j) = airev_2d(i, j) / (cuv(i, j) * cuv(i, j))  
           aivscu2gam_2d(i, j) = airvscu2_2d(i, j)**(-gamdi_grot)  
        END DO  
        airvscu2_2d(iip1, j) = airvscu2_2d(1, j)  
        aivscu2gam_2d(iip1, j) = aivscu2gam_2d(1, j)  
     END DO  
   
     DO j = 2, jjm  
        DO i = 1, iim  
           airuscv2_2d(i, j) = aireu_2d(i, j) / (cvu(i, j) * cvu(i, j))  
           aiuscv2gam_2d(i, j) = airuscv2_2d(i, j)**(-gamdi_grot)  
        END DO  
        airuscv2_2d(iip1, j) = airuscv2_2d(1, j)  
        aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j)  
     END DO  
375    
376      ! Calcul des aires aux pôles :      ! Calcul des aires aux pôles :
377    
# Line 534  contains Line 403  contains
403         constang_2d(iip1, j) = constang_2d(1, j)         constang_2d(iip1, j) = constang_2d(1, j)
404      END DO      END DO
405    
     call new_unit(unit)  
     open(unit, file="longitude_latitude.txt", status="replace", action="write")  
     write(unit, fmt=*) '"longitudes at V points (degrees)"', rlonv * 180. / pi  
     write(unit, fmt=*) '"latitudes at V points (degrees)"', rlatv * 180. / pi  
     write(unit, fmt=*) '"longitudes at U points (degrees)"', rlonu * 180. / pi  
     write(unit, fmt=*) '"latitudes at U points (degrees)"', rlatu * 180. / pi  
     close(unit)  
   
406    END SUBROUTINE inigeom    END SUBROUTINE inigeom
407    
408  end module comgeom  end module comgeom
Legend:
Removed from v.96  
changed lines
  Added in v.178
  ViewVC Help
Powered by ViewVC 1.1.21