/[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

revision 134 by guez, Wed Apr 29 15:47:56 2015 UTC 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      ! Fonction "f(y)" à dérivée tangente hyperbolique. Calcul des      ! Calcul des coefficients cu_2d, cv_2d, 1. / cu_2d**2, 1. /
115      ! coefficients cu_2d, cv_2d, 1. / cu_2d**2, 1. / cv_2d**2. Les      ! cv_2d**2. Les coefficients cu_2d et cv_2d permettent de passer
116      ! coefficients cu_2d et cv_2d permettent de passer des vitesses      ! des vitesses naturelles aux vitesses covariantes et
117      ! naturelles aux vitesses covariantes et contravariantes, ou      ! contravariantes, ou vice-versa.
     ! 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 fxhyp_m, only: fxhyp      use dynetat0_m, only: xprimp025, xprimm025, rlatu1, rlatu2, rlatu, rlatv, &
143      use fyhyp_m, only: fyhyp           yprimu1, yprimu2
     use jumble, only: new_unit  
     use nr_util, only: pi  
144      USE paramet_m, ONLY : iip1, jjp1      USE paramet_m, ONLY : iip1, jjp1
145    
146      ! Local:      ! Local:
147      INTEGER i, j, unit      INTEGER i, j
     REAL cvu(iip1, jjp1), cuv(iip1, jjm)  
148      REAL ai14, ai23, airez, un4rad2      REAL ai14, ai23, airez, un4rad2
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 yprimu(jjp1)  
152      REAL gamdi_gdiv, gamdi_grot, gamdi_h      REAL gamdi_gdiv, gamdi_grot, gamdi_h
     REAL xprimm025(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    
# Line 212  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    
     CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)  
     CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)  
   
     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 365  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 385  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 419  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 473  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.134  
changed lines
  Added in v.178
  ViewVC Help
Powered by ViewVC 1.1.21