trunk/dyn3d/tau2alpha.f

module tau2alpha_m

  USE paramet_m, ONLY : iip1, jjp1
  USE dimens_m, ONLY : jjm

  IMPLICIT NONE

  private iip1, jjp1, jjm

  REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)

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
    USE nr_util, ONLY : pi
    USE serre, ONLY : clat, clon, grossismx, grossismy

    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

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

    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
       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
       DO j = 1, jjm
          DO i = 1, iip1
             dxdyv(i, j) = 0.5 * (dxdys(i, j) + dxdys(i, j + 1))
          END DO
       END DO

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