/[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 161 by guez, Fri Jul 24 14:27:59 2015 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  
144      use nr_util, only: pi      use nr_util, only: pi
145      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  
146    
147      ! Local:      ! Local:
148      INTEGER i, j, itmax, itmay, iter, unit      INTEGER i, j
     REAL cvu(iip1, jjp1), cuv(iip1, jjm)  
149      REAL ai14, ai23, airez, un4rad2      REAL ai14, ai23, airez, un4rad2
     REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm  
150      REAL coslatm, coslatp, radclatm, radclatp      REAL coslatm, coslatp, radclatm, radclatp
151      REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m      REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
152      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)  
153      REAL gamdi_gdiv, gamdi_grot, gamdi_h      REAL gamdi_gdiv, gamdi_grot, gamdi_h
     REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)  
154      real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &      real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
155           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)  
156    
157      !------------------------------------------------------------------      !------------------------------------------------------------------
158    
159      PRINT *, 'Call sequence information: inigeom'      PRINT *, 'Call sequence information: inigeom'
160    
161      IF (nitergdiv/=2) THEN      IF (nitergdiv /= 2) THEN
162         gamdi_gdiv = coefdis / (real(nitergdiv)-2.)         gamdi_gdiv = coefdis / (nitergdiv - 2)
163      ELSE      ELSE
164         gamdi_gdiv = 0.         gamdi_gdiv = 0.
165      END IF      END IF
166      IF (nitergrot/=2) THEN  
167         gamdi_grot = coefdis / (real(nitergrot)-2.)      IF (nitergrot /= 2) THEN
168           gamdi_grot = coefdis / (nitergrot - 2)
169      ELSE      ELSE
170         gamdi_grot = 0.         gamdi_grot = 0.
171      END IF      END IF
172      IF (niterh/=2) THEN  
173         gamdi_h = coefdis / (real(niterh)-2.)      IF (niterh /= 2) THEN
174           gamdi_h = coefdis / (niterh - 2)
175      ELSE      ELSE
176         gamdi_h = 0.         gamdi_h = 0.
177      END IF      END IF
# Line 215  contains Line 180  contains
180      print *, "gamdi_grot = ", gamdi_grot      print *, "gamdi_grot = ", gamdi_grot
181      print *, "gamdi_h = ", gamdi_h      print *, "gamdi_h = ", gamdi_h
182    
     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.  
   
183      un4rad2 = 0.25 * rad * rad      un4rad2 = 0.25 * rad * rad
184    
185      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
# Line 426  contains Line 322  contains
322         unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)         unsairz_gam_2d(iip1, j) = unsairz_gam_2d(1, j)
323      END DO      END DO
324    
325      ! Calcul des élongations cu_2d, cv_2d, cvu      ! Calcul des élongations cu_2d, cv_2d
326    
327      DO j = 1, jjm      DO j = 1, jjm
328         DO i = 1, iim         DO i = 1, iim
329            cv_2d(i, j) = 0.5 * &            cv_2d(i, j) = 0.5 * &
330                 (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))  
331            unscv2_2d(i, j) = 1. / cv_2d(i, j)**2            unscv2_2d(i, j) = 1. / cv_2d(i, j)**2
332         END DO         END DO
333         DO i = 1, iim         DO i = 1, iim
# Line 446  contains Line 338  contains
338            cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)            cvscuvgam_2d(i, j) = cvsurcuv_2d(i, j)**(-gamdi_grot)
339         END DO         END DO
340         cv_2d(iip1, j) = cv_2d(1, j)         cv_2d(iip1, j) = cv_2d(1, j)
        cvu(iip1, j) = cvu(1, j)  
341         unscv2_2d(iip1, j) = unscv2_2d(1, j)         unscv2_2d(iip1, j) = unscv2_2d(1, j)
        cuv(iip1, j) = cuv(1, j)  
342         cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)         cuvsurcv_2d(iip1, j) = cuvsurcv_2d(1, j)
343         cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)         cvsurcuv_2d(iip1, j) = cvsurcuv_2d(1, j)
344         cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)         cuvscvgam1_2d(iip1, j) = cuvscvgam1_2d(1, j)
# Line 480  contains Line 370  contains
370    
371      cu_2d(:, 1) = 0.      cu_2d(:, 1) = 0.
372      unscu2_2d(:, 1) = 0.      unscu2_2d(:, 1) = 0.
     cvu(:, 1) = 0.  
373    
374      cu_2d(:, jjp1) = 0.      cu_2d(:, jjp1) = 0.
375      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  
376    
377      ! Calcul des aires aux pôles :      ! Calcul des aires aux pôles :
378    
# Line 534  contains Line 404  contains
404         constang_2d(iip1, j) = constang_2d(1, j)         constang_2d(iip1, j) = constang_2d(1, j)
405      END DO      END DO
406    
     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)  
   
407    END SUBROUTINE inigeom    END SUBROUTINE inigeom
408    
409  end module comgeom  end module comgeom
Legend:
Removed from v.96  
changed lines
  Added in v.161
  ViewVC Help
Powered by ViewVC 1.1.21