trunk/dyn3d/tau2alpha.f

module tau2alpha_m

  IMPLICIT NONE

contains

  SUBROUTINE tau2alpha(type, factt, taumin, taumax, alpha)

    USE comgeom, ONLY: cu_2d, cv_2d, rlatu, rlatv
    use conf_guide_m, only: lat_min_guide, lat_max_guide
    USE dimens_m, ONLY: iim, jjm
    USE nr_util, ONLY: pi
    USE paramet_m, ONLY: iip1, jjp1
    USE serre, ONLY: clat, clon, grossismx, grossismy
    use writefield_m, only: writefield

    INTEGER, intent(in):: type
    REAL, intent(in):: factt, taumin, taumax
    real, intent(out):: alpha(:, :)

    ! Local:
    REAL dxdy
    REAL dxdy_min, dxdy_max
    REAL alphamin, alphamax, xi
    REAL, SAVE:: gamma
    INTEGER i, j, ilon, ilat
    LOGICAL:: first = .TRUE.
    REAL zdx(iip1, jjp1), zdy(iip1, jjp1)
    REAL zlat
    REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)

    !------------------------------------------------------------

    IF (first) THEN
       DO j = 2, jjm
          DO i = 2, iip1
             zdx(i, j) = 0.5 * (cu_2d(i - 1, j) + cu_2d(i, j)) / cos(rlatu(j))
          END DO
          zdx(1, j) = zdx(iip1, j)
       END DO
       DO j = 2, jjm
          DO i = 1, iip1
             zdy(i, j) = 0.5 * (cv_2d(i, j - 1) + cv_2d(i, j))
          END DO
       END DO
       DO i = 1, iip1
          zdx(i, 1) = zdx(i, 2)
          zdx(i, jjp1) = zdx(i, jjm)
          zdy(i, 1) = zdy(i, 2)
          zdy(i, jjp1) = zdy(i, jjm)
       END DO

       DO j = 1, jjp1
          DO i = 1, iip1
             dxdys(i, j) = sqrt(zdx(i, j)**2 + zdy(i, j)**2)
          END DO
       END DO
       CALL writefield("dxdys", dxdys)

       if (type == 2) then
          DO j = 1, jjp1
             DO i = 1, iim
                dxdyu(i, j) = 0.5 * (dxdys(i, j) + dxdys(i + 1, j))
             END DO
             dxdyu(iip1, j) = dxdyu(1, j)
          END DO
       elseif (type == 3) then
          DO j = 1, jjm
             DO i = 1, iip1
                dxdyv(i, j) = 0.5 * (dxdys(i, j) + dxdys(i, j + 1))
             END DO
          END DO
       end if

       ! coordonnees du centre du zoom
       CALL coordij(clon, clat, ilon, ilat)
       ! aire de la maille au centre du zoom
       dxdy_min = dxdys(ilon, ilat)

       ! dxdy maximal de la maille :
       dxdy_max = 0.
       DO j = 1, jjp1
          DO i = 1, iip1
             dxdy_max = max(dxdy_max, dxdys(i, j))
          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.'
       ELSE
          gamma = (dxdy_max - 2. * dxdy_min) / (dxdy_max - dxdy_min)
          PRINT *, 'gamma=', gamma
          IF (gamma<1.E-5) THEN
             PRINT *, 'gamma =', gamma, '<1e-5'
             STOP
          END IF
          PRINT *, 'gamma=', gamma
          gamma = log(0.5) / log(gamma)
       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
    END DO

  END SUBROUTINE tau2alpha

end module tau2alpha_m
1	module tau2alpha_m
2
3	IMPLICIT NONE
4
5	contains
6
7	SUBROUTINE tau2alpha(type, factt, taumin, taumax, alpha)
8
9	USE comgeom, ONLY: cu_2d, cv_2d, rlatu, rlatv
10	use conf_guide_m, only: lat_min_guide, lat_max_guide
11	USE dimens_m, ONLY: iim, jjm
12	USE nr_util, ONLY: pi
13	USE paramet_m, ONLY: iip1, jjp1
14	USE serre, ONLY: clat, clon, grossismx, grossismy
15	use writefield_m, only: writefield
16
17	INTEGER, intent(in):: type
18	REAL, intent(in):: factt, taumin, taumax
19	real, intent(out):: alpha(:, :)
20
21	! Local:
22	REAL dxdy
23	REAL dxdy_min, dxdy_max
24	REAL alphamin, alphamax, xi
25	REAL, SAVE:: gamma
26	INTEGER i, j, ilon, ilat
27	LOGICAL:: first = .TRUE.
28	REAL zdx(iip1, jjp1), zdy(iip1, jjp1)
29	REAL zlat
30	REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
31
32	!------------------------------------------------------------
33
34	IF (first) THEN
35	DO j = 2, jjm
36	DO i = 2, iip1
37	zdx(i, j) = 0.5 * (cu_2d(i - 1, j) + cu_2d(i, j)) / cos(rlatu(j))
38	END DO
39	zdx(1, j) = zdx(iip1, j)
40	END DO
41	DO j = 2, jjm
42	DO i = 1, iip1
43	zdy(i, j) = 0.5 * (cv_2d(i, j - 1) + cv_2d(i, j))
44	END DO
45	END DO
46	DO i = 1, iip1
47	zdx(i, 1) = zdx(i, 2)
48	zdx(i, jjp1) = zdx(i, jjm)
49	zdy(i, 1) = zdy(i, 2)
50	zdy(i, jjp1) = zdy(i, jjm)
51	END DO
52
53	DO j = 1, jjp1
54	DO i = 1, iip1
55	dxdys(i, j) = sqrt(zdx(i, j)2 + zdy(i, j)2)
56	END DO
57	END DO
58	CALL writefield("dxdys", dxdys)
59
60	if (type == 2) then
61	DO j = 1, jjp1
62	DO i = 1, iim
63	dxdyu(i, j) = 0.5 * (dxdys(i, j) + dxdys(i + 1, j))
64	END DO
65	dxdyu(iip1, j) = dxdyu(1, j)
66	END DO
67	elseif (type == 3) then
68	DO j = 1, jjm
69	DO i = 1, iip1
70	dxdyv(i, j) = 0.5 * (dxdys(i, j) + dxdys(i, j + 1))
71	END DO
72	END DO
73	end if
74
75	! coordonnees du centre du zoom
76	CALL coordij(clon, clat, ilon, ilat)
77	! aire de la maille au centre du zoom
78	dxdy_min = dxdys(ilon, ilat)
79
80	! dxdy maximal de la maille :
81	dxdy_max = 0.
82	DO j = 1, jjp1
83	DO i = 1, iip1
84	dxdy_max = max(dxdy_max, dxdys(i, j))
85	END DO
86	END DO
87
88	IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN
89	PRINT *, 'Attention : modèle peu zoomé.'
90	PRINT *, 'On prend une constante de guidage constante.'
91	ELSE
92	gamma = (dxdy_max - 2. * dxdy_min) / (dxdy_max - dxdy_min)
93	PRINT *, 'gamma=', gamma
94	IF (gamma<1.E-5) THEN
95	PRINT *, 'gamma =', gamma, '<1e-5'
96	STOP
97	END IF
98	PRINT *, 'gamma=', gamma
99	gamma = log(0.5) / log(gamma)
100	END IF
101	first = .false.
102	END IF
103
104	alphamin = factt / taumax
105	alphamax = factt / taumin
106
107	DO j = 1, size(alpha, 2)
108	DO i = 1, size(alpha, 1)
109	IF (type==1) THEN
110	dxdy = dxdys(i, j)
111	zlat = rlatu(j) * 180. / pi
112	ELSE IF (type==2) THEN
113	dxdy = dxdyu(i, j)
114	zlat = rlatu(j) * 180. / pi
115	ELSE IF (type==3) THEN
116	dxdy = dxdyv(i, j)
117	zlat = rlatv(j) * 180. / pi
118	END IF
119	IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN
120	! grille regulière
121	alpha(i, j) = alphamin
122	ELSE
123	xi = ((dxdy_max - dxdy) / (dxdy_max - dxdy_min))**gamma
124	xi = min(xi, 1.)
125	IF (lat_min_guide <= zlat .AND. zlat <= lat_max_guide) THEN
126	alpha(i, j) = xi * alphamin + (1. - xi) * alphamax
127	ELSE
128	alpha(i, j) = 0.
129	END IF
130	END IF
131	END DO
132	END DO
133
134	END SUBROUTINE tau2alpha
135
136	end module tau2alpha_m