/[lmdze]/trunk/Sources/dyn3d/Guide/tau2alpha.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/Guide/tau2alpha.f

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

revision 114 by guez, Fri Sep 19 11:41:35 2014 UTC revision 115 by guez, Fri Sep 19 17:36:20 2014 UTC
# Line 4  module tau2alpha_m Line 4  module tau2alpha_m
4    
5  contains  contains
6    
7    SUBROUTINE tau2alpha(type, factt, taumin, taumax, alpha)    SUBROUTINE tau2alpha(dxdy, rlat, taumin, taumax, alpha)
8    
9      USE comgeom, ONLY: cu_2d, cv_2d, rlatu, rlatv      use conf_guide_m, only: lat_min_guide, lat_max_guide, factt
10      use conf_guide_m, only: lat_min_guide, lat_max_guide      use init_tau2alpha_m, only: dxdy_min, dxdy_max, gamma
11      USE dimens_m, ONLY: iim, jjm      USE nr_util, ONLY: pi, assert_eq
12      USE nr_util, ONLY: pi  
13      USE paramet_m, ONLY: iip1, jjp1      REAL, intent(in):: dxdy(:, :) ! (n_lon, n_lat)
14      USE serre, ONLY: clat, clon, grossismx, grossismy      REAL, intent(in):: rlat(:) ! (n_lat)
15      use writefield_m, only: writefield      REAL, intent(in):: taumin, taumax
16        real, intent(out):: alpha(:, :) ! (n_lon, n_lat)
     INTEGER, intent(in):: type  
     REAL, intent(in):: factt, taumin, taumax  
     real, intent(out):: alpha(:, :)  
17    
18      ! Local:      ! Local:
     REAL dxdy  
     REAL, save:: dxdy_min, dxdy_max  
19      REAL alphamin, alphamax, xi      REAL alphamin, alphamax, xi
20      REAL, SAVE:: gamma      INTEGER i, j, n_lon, n_lat
     INTEGER i, j, ilon, ilat  
     LOGICAL:: first = .TRUE.  
     REAL dx(iip1, jjp1), dy(iip1, jjp1)  
21      REAL zlat      REAL zlat
     REAL, save:: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)  
22    
23      !------------------------------------------------------------      !------------------------------------------------------------
24    
25      IF (first) THEN      PRINT *, 'Call sequence information: tau2alpha'
        DO j = 2, jjm  
           DO i = 2, iip1  
              dx(i, j) = 0.5 * (cu_2d(i - 1, j) + cu_2d(i, j)) / cos(rlatu(j))  
           END DO  
           dx(1, j) = dx(iip1, j)  
        END DO  
        DO j = 2, jjm  
           DO i = 1, iip1  
              dy(i, j) = 0.5 * (cv_2d(i, j - 1) + cv_2d(i, j))  
           END DO  
        END DO  
        DO i = 1, iip1  
           dx(i, 1) = dx(i, 2)  
           dx(i, jjp1) = dx(i, jjm)  
           dy(i, 1) = dy(i, 2)  
           dy(i, jjp1) = dy(i, jjm)  
        END DO  
   
        DO j = 1, jjp1  
           DO i = 1, iip1  
              dxdys(i, j) = sqrt(dx(i, j)**2 + dy(i, j)**2)  
           END DO  
        END DO  
        CALL writefield("dxdys", dxdys)  
26    
27         DO j = 1, jjp1      n_lon = assert_eq(size(alpha, 1), size(dxdy, 1), "tau2alpha n_lon")
28            DO i = 1, iim      n_lat = assert_eq(size(alpha, 2), size(dxdy, 2), size(rlat), &
29               dxdyu(i, j) = 0.5 * (dxdys(i, j) + dxdys(i + 1, j))           "tau2alpha n_lat")
           END DO  
           dxdyu(iip1, j) = dxdyu(1, j)  
        END DO  
30    
31         DO j = 1, jjm      alphamin = factt / taumax
32            DO i = 1, iip1      alphamax = factt / taumin
              dxdyv(i, j) = 0.5 * (dxdys(i, j) + dxdys(i, j + 1))  
           END DO  
        END DO  
33    
34         ! coordonnees du centre du zoom      DO j = 1, n_lat
35         CALL coordij(clon, clat, ilon, ilat)         zlat = rlat(j) * 180. / pi
36         ! aire de la maille au centre du zoom         IF (lat_min_guide <= zlat .AND. zlat <= lat_max_guide) THEN
37         dxdy_min = dxdys(ilon, ilat)            DO i = 1, n_lon
38                 xi = min(((dxdy_max - dxdy(i, j)) &
39         ! dxdy maximal de la maille :                    / (dxdy_max - dxdy_min))**gamma, 1.)
40         dxdy_max = 0.               alpha(i, j) = xi * alphamin + (1. - xi) * alphamax
        DO j = 1, jjp1  
           DO i = 1, iip1  
              dxdy_max = max(dxdy_max, dxdys(i, j))  
41            END DO            END DO
        END DO  
   
        IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN  
           PRINT *, 'Attention : modèle peu zoomé.'  
           PRINT *, 'On prend une constante de guidage constante.'  
42         ELSE         ELSE
43            gamma = (dxdy_max - 2. * dxdy_min) / (dxdy_max - dxdy_min)            alpha(:, j) = 0.
           IF (gamma < 1E-5) THEN  
              PRINT *, '(dxdy_max - 2. * dxdy_min) / (dxdy_max - dxdy_min) ' &  
                   // '< 1e-5'  
              STOP 1  
           END IF  
           gamma = log(0.5) / log(gamma)  
           PRINT *, 'gamma=', gamma  
44         END IF         END IF
        first = .false.  
     END IF  
   
     alphamin = factt / taumax  
     alphamax = factt / taumin  
   
     DO j = 1, size(alpha, 2)  
        DO i = 1, size(alpha, 1)  
           IF (type==1) THEN  
              dxdy = dxdys(i, j)  
              zlat = rlatu(j) * 180. / pi  
           ELSE IF (type==2) THEN  
              dxdy = dxdyu(i, j)  
              zlat = rlatu(j) * 180. / pi  
           ELSE IF (type==3) THEN  
              dxdy = dxdyv(i, j)  
              zlat = rlatv(j) * 180. / pi  
           END IF  
           IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN  
              ! grille regulière  
              alpha(i, j) = alphamin  
           ELSE  
              xi = ((dxdy_max - dxdy) / (dxdy_max - dxdy_min))**gamma  
              xi = min(xi, 1.)  
              IF (lat_min_guide <= zlat .AND. zlat <= lat_max_guide) THEN  
                 alpha(i, j) = xi * alphamin + (1. - xi) * alphamax  
              ELSE  
                 alpha(i, j) = 0.  
              END IF  
           END IF  
        END DO  
45      END DO      END DO
46    
47    END SUBROUTINE tau2alpha    END SUBROUTINE tau2alpha
Legend:
Removed from v.114  
changed lines
  Added in v.115
  ViewVC Help
Powered by ViewVC 1.1.21