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

Diff of /trunk/dyn3d/comgeom.f

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

trunk/dyn3d/comgeom.f revision 97 by guez, Fri Apr 25 14:58:31 2014 UTC trunk/Sources/dyn3d/comgeom.f revision 139 by guez, Tue May 26 17:46:03 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 à      ! Fonction "f(y)" à dérivée tangente hyperbolique. Calcul des
115      ! dérivée tangente hyperbolique. Calcul des coefficients cu_2d,      ! coefficients cu_2d, cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les
116      ! cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les coefficients cu_2d et cv_2d      ! coefficients cu_2d et cv_2d permettent de passer des vitesses
117      ! permettent de passer des vitesses naturelles aux vitesses      ! naturelles aux vitesses covariantes et contravariantes, ou
118      ! covariantes et contravariantes, ou vice-versa.      ! vice-versa.
119    
120      ! On a :      ! On a :
121      ! 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 152  contains Line 139  contains
139      ! 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
140      ! adressage plus facile en ij.      ! adressage plus facile en ij.
141    
142      ! xprimu et xprimv sont respectivement les valeurs de dx / dX aux      ! cvu et cv_2d sont respectivement les valeurs de
143      ! points u et v. yprimu et yprimv sont respectivement les valeurs      ! cv_2d aux points u et v.
     ! 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.  
144    
145      ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d      ! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d
146      ! aux points u, v, scalaires, et z. Cf. "inigeom.txt".      ! aux points u, v, scalaires, et z. Cf. "inigeom.txt".
147    
148      USE comconst, ONLY : g, omeg, rad      USE comconst, ONLY : g, omeg, rad
149      USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh      USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh
150      use conf_gcm_m, ONLY : fxyhypb, ysinus      use dynetat0_m, only: xprimp025, xprimm025, rlatu1, rlatu2, rlatu, rlatv, &
151      use fxy_m, only: fxy           yprimu1, yprimu2, rlonu, rlonv
     use fxyhyper_m, only: fxyhyper  
     use fxysinus_m, only: fxysinus  
152      use jumble, only: new_unit      use jumble, only: new_unit
153      use nr_util, only: pi      use nr_util, only: pi
154      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  
155    
156      ! Local:      ! Local:
157      INTEGER i, j, itmax, itmay, iter, unit      INTEGER i, j, unit
158      REAL cvu(iip1, jjp1), cuv(iip1, jjm)      REAL cvu(iip1, jjp1), cuv(iip1, jjm)
159      REAL ai14, ai23, airez, un4rad2      REAL ai14, ai23, airez, un4rad2
     REAL eps, x1, xo1, f, df, xdm, y1, yo1, ydm  
160      REAL coslatm, coslatp, radclatm, radclatp      REAL coslatm, coslatp, radclatm, radclatp
161      REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m      REAL, dimension(iip1, jjp1):: cuij1, cuij2, cuij3, cuij4 ! in m
162      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)  
163      REAL gamdi_gdiv, gamdi_grot, gamdi_h      REAL gamdi_gdiv, gamdi_grot, gamdi_h
     REAL rlonm025(iip1), xprimm025(iip1), rlonp025(iip1), xprimp025(iip1)  
164      real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &      real, dimension(iim + 1, jjm + 1):: aireij1_2d, aireij2_2d, aireij3_2d, &
165           aireij4_2d ! in m2           aireij4_2d ! in m2
166      real airuscv2_2d(iim + 1, jjm)      real airuscv2_2d(iim + 1, jjm)
# Line 196  contains Line 171  contains
171    
172      PRINT *, 'Call sequence information: inigeom'      PRINT *, 'Call sequence information: inigeom'
173    
174      IF (nitergdiv/=2) THEN      IF (nitergdiv /= 2) THEN
175         gamdi_gdiv = coefdis / (real(nitergdiv)-2.)         gamdi_gdiv = coefdis / (nitergdiv - 2)
176      ELSE      ELSE
177         gamdi_gdiv = 0.         gamdi_gdiv = 0.
178      END IF      END IF
179      IF (nitergrot/=2) THEN  
180         gamdi_grot = coefdis / (real(nitergrot)-2.)      IF (nitergrot /= 2) THEN
181           gamdi_grot = coefdis / (nitergrot - 2)
182      ELSE      ELSE
183         gamdi_grot = 0.         gamdi_grot = 0.
184      END IF      END IF
185      IF (niterh/=2) THEN  
186         gamdi_h = coefdis / (real(niterh)-2.)      IF (niterh /= 2) THEN
187           gamdi_h = coefdis / (niterh - 2)
188      ELSE      ELSE
189         gamdi_h = 0.         gamdi_h = 0.
190      END IF      END IF
# Line 216  contains Line 193  contains
193      print *, "gamdi_grot = ", gamdi_grot      print *, "gamdi_grot = ", gamdi_grot
194      print *, "gamdi_h = ", gamdi_h      print *, "gamdi_h = ", gamdi_h
195    
     IF (fxyhypb) THEN  
        ! 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)  
     ELSE  
        IF (ysinus) THEN  
           print *, 'inigeom: Y = sin(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  
     END IF  
   
     rlatu(1) = pi / 2.  
     rlatu(jjp1) = -rlatu(1)  
   
     ! Calcul aux pôles  
   
     yprimu(1) = 0.  
     yprimu(jjp1) = 0.  
   
196      un4rad2 = 0.25 * rad * rad      un4rad2 = 0.25 * rad * rad
197    
198      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires      ! Cf. "inigeom.txt". Calcul des quatre aires élémentaires
Legend:
Removed from v.97  
changed lines
  Added in v.139
  ViewVC Help
Powered by ViewVC 1.1.21