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 dump2d_m, only: dump2d
    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

       CALL dump2d(iip1, jjp1, dxdys, 'DX2DY2 SCAL ')
       CALL dump2d(iip1, jjp1, dxdyu, 'DX2DY2 U ')
       CALL dump2d(iip1, jjp1, dxdyv, 'DX2DY2 v ')

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