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	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	102	use dump2d_m, only: dump2d
20	guez	39	USE nr_util, ONLY : pi
21	guez	37	USE serre, ONLY : clat, clon, grossismx, grossismy
22
23	guez	103	INTEGER, intent(in):: type
24	guez	44	REAL, intent(in):: factt, taumin, taumax
25	guez	103	real, intent(out):: alpha(:, :)
26
27			! Local:
28			REAL dxdy
29	guez	37	REAL dxdy_min, dxdy_max
30	guez	103	REAL alphamin, alphamax, xi
31			REAL, SAVE:: gamma
32	guez	37	INTEGER i, j, ilon, ilat
33	guez	103	LOGICAL:: first = .TRUE.
34	guez	37	REAL zdx(iip1, jjp1), zdy(iip1, jjp1)
35			REAL zlat
36
37	guez	44	!------------------------------------------------------------
38
39	guez	37	IF (first) THEN
40			DO j = 2, jjm
41			DO i = 2, iip1
42	guez	103	zdx(i, j) = 0.5 * (cu_2d(i - 1, j) + cu_2d(i, j)) / cos(rlatu(j))
43	guez	37	END DO
44			zdx(1, j) = zdx(iip1, j)
45			END DO
46			DO j = 2, jjm
47			DO i = 1, iip1
48	guez	103	zdy(i, j) = 0.5 * (cv_2d(i, j - 1) + cv_2d(i, j))
49	guez	37	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	guez	103	dxdys(i, j) = sqrt(zdx(i, j)2 + zdy(i, j)2)
60	guez	37	END DO
61			END DO
62			DO j = 1, jjp1
63			DO i = 1, iim
64	guez	103	dxdyu(i, j) = 0.5 * (dxdys(i, j) + dxdys(i + 1, j))
65	guez	37	END DO
66			dxdyu(iip1, j) = dxdyu(1, j)
67			END DO
68			DO j = 1, jjm
69			DO i = 1, iip1
70	guez	103	dxdyv(i, j) = 0.5 * (dxdys(i, j) + dxdys(i, j + 1))
71	guez	37	END DO
72			END DO
73
74	guez	103	CALL dump2d(iip1, jjp1, dxdys, 'DX2DY2 SCAL ')
75			CALL dump2d(iip1, jjp1, dxdyu, 'DX2DY2 U ')
76			CALL dump2d(iip1, jjp1, dxdyv, 'DX2DY2 v ')
77	guez	37
78	guez	103	! coordonnees du centre du zoom
79	guez	37	CALL coordij(clon, clat, ilon, ilat)
80	guez	103	! aire de la maille au centre du zoom
81	guez	37	dxdy_min = dxdys(ilon, ilat)
82	guez	103
83			! dxdy maximal de la maille :
84	guez	37	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	guez	103	IF (abs(grossismx - 1.)<0.1 .OR. abs(grossismy - 1.)<0.1) THEN
92	guez	37	PRINT *, 'ATTENTION modele peu zoome'
93			PRINT *, 'ATTENTION on prend une constante de guidage cste'
94			gamma = 0.
95			ELSE
96	guez	103	gamma = (dxdy_max - 2. * dxdy_min) / (dxdy_max - dxdy_min)
97	guez	37	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	guez	103	gamma = log(0.5) / log(gamma)
104	guez	37	END IF
105	guez	103	first = .false.
106	guez	37	END IF
107
108	guez	103	alphamin = factt / taumax
109			alphamax = factt / taumin
110	guez	37
111	guez	103	DO j = 1, size(alpha, 2)
112			DO i = 1, size(alpha, 1)
113	guez	37	IF (type==1) THEN
114	guez	103	dxdy = dxdys(i, j)
115			zlat = rlatu(j) * 180. / pi
116	guez	37	ELSE IF (type==2) THEN
117	guez	103	dxdy = dxdyu(i, j)
118			zlat = rlatu(j) * 180. / pi
119	guez	37	ELSE IF (type==3) THEN
120	guez	103	dxdy = dxdyv(i, j)
121			zlat = rlatv(j) * 180. / pi
122	guez	37	END IF
123	guez	103	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	guez	37	alpha(i, j) = alphamin
126			ELSE
127	guez	103	xi = ((dxdy_max - dxdy) / (dxdy_max - dxdy_min))**gamma
128	guez	37	xi = min(xi, 1.)
129	guez	103	IF (lat_min_guide <= zlat .AND. zlat <= lat_max_guide) THEN
130			alpha(i, j) = xi * alphamin + (1. - xi) * alphamax
131	guez	37	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