dyn3d/Guide/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, save:: dxdy_min, dxdy_max
    REAL alphamin, alphamax, xi
    REAL, SAVE:: gamma
    INTEGER i, j, ilon, ilat
    LOGICAL:: first = .TRUE.
    REAL dx(iip1, jjp1), dy(iip1, jjp1)
    REAL zlat
    REAL, save:: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)

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

    IF (first) THEN
       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)

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