dyn3d/Guide/tau2alpha.f

module tau2alpha_m

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

    IMPLICIT NONE

    private iip1, jjp1, jjm

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

contains

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

    USE dimens_m, ONLY : iim
    USE nr_util, ONLY : pi
    USE comgeom, ONLY : cu_2d, cv_2d, rlatu, rlatv
    USE serre, ONLY : clat, clon, grossismx, grossismy

    !   arguments :
    INTEGER type
    INTEGER pim, pjm
    REAL, intent(in):: factt, taumin, taumax
    REAL dxdy_, alpha(pim, pjm)
    REAL dxdy_min, dxdy_max

    !  local :
    REAL alphamin, alphamax, gamma, xi
    SAVE gamma
    INTEGER i, j, ilon, ilat

    LOGICAL first
    SAVE first
    DATA 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)*zdx(i, j)+zdy(i, j)*zdy(i, j))
          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+1, j))
          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 maximale 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
    END IF

    alphamin = factt/taumax
    alphamax = factt/taumin

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