dyn3d/Dissipation/inidissip.f

module inidissip_m

  use dimens_m, only: llm

  IMPLICIT NONE

  private llm

  REAL dtdiss
  integer idissip ! periode de la dissipation (en pas)
  real tetaudiv(llm),tetaurot(llm),tetah(llm)   
  real cdivu, crot, cdivh

contains

  SUBROUTINE inidissip

    ! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06
    ! Initialisation de la dissipation horizontale                        

    USE comconst, ONLY : dtvr
    use comdissnew, only: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
         tetagrot, tetatemp
    USE comvert, ONLY : preff, presnivs
    USE conf_gcm_m, ONLY : iperiod
    USE dimens_m, ONLY : jjm, llm
    USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
    use new_unit_m, only: new_unit
    use filtreg_m, only: filtreg

    ! Variables local to the procedure:
    REAL zvert(llm), max_zvert
    REAL zh(ip1jmp1), zu(ip1jmp1), zv(ip1jm), deltap(ip1jmp1, llm)
    REAL ullm, vllm, umin, vmin, zhmin, zhmax
    REAL zllm, z1llm
    INTEGER l, ij, idum, ii, unit
    REAL tetamin ! in s
    REAL ran1

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

    PRINT *, 'Call sequence information: inidissip'

    !   calcul des valeurs propres des operateurs par methode iterrative:   

    crot = -1.
    cdivu = -1.
    cdivh = -1.

    !   calcul de la valeur propre de divgrad:                              

    idum = 0
    DO l = 1, llm
       DO ij = 1, ip1jmp1
          deltap(ij, l) = 1.
       END DO
    END DO

    idum = -1
    zh(1) = ran1(idum) - .5
    idum = 0
    DO ij = 2, ip1jmp1
       zh(ij) = ran1(idum) - .5
    END DO

    CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)

    CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
    IF (zhmin >= zhmax) THEN
       PRINT *, 'zhmin zhmax', zhmin, zhmax
       print *, 'Problème générateur aléatoire dans inidissip'
       STOP 1
    END IF

    zllm = abs(zhmax)
    DO l = 1, 50
       IF (lstardis) THEN
          CALL divgrad2(1, zh, deltap, niterh, zh)
       ELSE
          CALL divgrad(1, zh, niterh, zh)
       END IF

       CALL minmax(iip1*jjp1, zh, zhmin, zhmax)

       zllm = abs(zhmax)
       z1llm = 1./zllm
       DO ij = 1, ip1jmp1
          zh(ij) = zh(ij)*z1llm
       END DO
    END DO

    IF (lstardis) THEN
       cdivh = 1./zllm
    ELSE
       cdivh = zllm**(-1./niterh)
    END IF

    !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)     

    PRINT *, 'calcul des valeurs propres'

    DO  ii = 1, 2
       DO ij = 1, ip1jmp1
          zu(ij) = ran1(idum) - .5
       END DO
       CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
       DO ij = 1, ip1jm
          zv(ij) = ran1(idum) - .5
       END DO
       CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)

       CALL minmax(iip1*jjp1, zu, umin, ullm)
       CALL minmax(iip1*jjm, zv, vmin, vllm)

       ullm = abs(ullm)
       vllm = abs(vllm)

       DO  l = 1, 50
          IF (ii==1) THEN
             IF (lstardis) THEN
                CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
             ELSE
                CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
             END IF
          ELSE
             IF (lstardis) THEN
                CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
             ELSE
                CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
             END IF
          END IF

          CALL minmax(iip1*jjp1, zu, umin, ullm)
          CALL minmax(iip1*jjm, zv, vmin, vllm)

          ullm = abs(ullm)
          vllm = abs(vllm)

          zllm = max(ullm, vllm)
          z1llm = 1./zllm
          DO ij = 1, ip1jmp1
             zu(ij) = zu(ij)*z1llm
          END DO
          DO ij = 1, ip1jm
             zv(ij) = zv(ij)*z1llm
          END DO
       end DO

       IF (ii==1) THEN
          IF (lstardis) THEN
             cdivu = 1./zllm
          ELSE
             cdivu = zllm**(-1./nitergdiv)
          END IF
       ELSE
          IF (lstardis) THEN
             crot = 1./zllm
          ELSE
             crot = zllm**(-1./nitergrot)
          END IF
       END IF
    END DO

    PRINT *, 'cdivu = ', cdivu
    PRINT *, 'crot = ', crot
    PRINT *, 'cdivh = ', cdivh

    ! Variation verticale du coefficient de dissipation :
    zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
    ! (between 1 and 2)

    tetaudiv = zvert / tetagdiv
    tetaurot = zvert / tetagrot
    tetah = zvert / tetatemp
    call new_unit(unit)
    open(unit, file="inidissip.csv", status="replace", action="write")
    write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
    do l = 1, llm
       write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
    end do
    close(unit)
    print *, 'Created file "inidissip.csv".'

    max_zvert = maxval(zvert)
    tetamin = min(1E6, tetagdiv / max_zvert, tetagrot / max_zvert, &
         tetatemp / max_zvert)
    PRINT *, 'tetamin = ', tetamin
    idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
    PRINT *, 'idissip = ', idissip
    dtdiss = idissip * dtvr
    PRINT *, 'dtdiss = ', dtdiss

  END SUBROUTINE inidissip

end module inidissip_m
1	guez	26	module inidissip_m
2	guez	3
3	guez	26	use dimens_m, only: llm
4	guez	3
5	guez	26	IMPLICIT NONE
6	guez	3
7	guez	26	private llm
8	guez	3
9	guez	26	REAL dtdiss
10			integer idissip ! periode de la dissipation (en pas)
11			real tetaudiv(llm),tetaurot(llm),tetah(llm)
12			real cdivu, crot, cdivh
13	guez	3
14	guez	26	contains
15	guez	3
16	guez	27	SUBROUTINE inidissip
17	guez	3
18	guez	26	! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06
19			! Initialisation de la dissipation horizontale
20	guez	3
21	guez	26	USE comconst, ONLY : dtvr
22	guez	27	use comdissnew, only: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
23			tetagrot, tetatemp
24	guez	26	USE comvert, ONLY : preff, presnivs
25			USE conf_gcm_m, ONLY : iperiod
26			USE dimens_m, ONLY : jjm, llm
27			USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1
28			use new_unit_m, only: new_unit
29	guez	27	use filtreg_m, only: filtreg
30	guez	3
31	guez	26	! Variables local to the procedure:
32	guez	27	REAL zvert(llm), max_zvert
33	guez	26	REAL zh(ip1jmp1), zu(ip1jmp1), zv(ip1jm), deltap(ip1jmp1, llm)
34			REAL ullm, vllm, umin, vmin, zhmin, zhmax
35			REAL zllm, z1llm
36			INTEGER l, ij, idum, ii, unit
37	guez	27	REAL tetamin ! in s
38	guez	26	REAL ran1
39	guez	3
40	guez	26	!-----------------------------------------------------------------------
41	guez	3
42	guez	26	PRINT *, 'Call sequence information: inidissip'
43	guez	3
44	guez	26	! calcul des valeurs propres des operateurs par methode iterrative:
45	guez	3
46	guez	26	crot = -1.
47			cdivu = -1.
48			cdivh = -1.
49	guez	3
50	guez	26	! calcul de la valeur propre de divgrad:
51	guez	3
52	guez	26	idum = 0
53			DO l = 1, llm
54			DO ij = 1, ip1jmp1
55			deltap(ij, l) = 1.
56			END DO
57			END DO
58	guez	3
59	guez	26	idum = -1
60			zh(1) = ran1(idum) - .5
61			idum = 0
62			DO ij = 2, ip1jmp1
63			zh(ij) = ran1(idum) - .5
64			END DO
65	guez	3
66	guez	26	CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
67	guez	3
68	guez	26	CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
69	guez	39	IF (zhmin >= zhmax) THEN
70			PRINT *, 'zhmin zhmax', zhmin, zhmax
71			print *, 'Problème générateur aléatoire dans inidissip'
72			STOP 1
73	guez	26	END IF
74	guez	3
75	guez	26	zllm = abs(zhmax)
76			DO l = 1, 50
77			IF (lstardis) THEN
78			CALL divgrad2(1, zh, deltap, niterh, zh)
79			ELSE
80			CALL divgrad(1, zh, niterh, zh)
81			END IF
82	guez	3
83	guez	26	CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
84	guez	3
85	guez	26	zllm = abs(zhmax)
86			z1llm = 1./zllm
87			DO ij = 1, ip1jmp1
88			zh(ij) = zh(ij)*z1llm
89			END DO
90			END DO
91	guez	3
92	guez	26	IF (lstardis) THEN
93			cdivh = 1./zllm
94			ELSE
95			cdivh = zllm**(-1./niterh)
96			END IF
97	guez	3
98	guez	26	! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2)
99	guez	3
100	guez	26	PRINT *, 'calcul des valeurs propres'
101	guez	3
102	guez	26	DO ii = 1, 2
103			DO ij = 1, ip1jmp1
104			zu(ij) = ran1(idum) - .5
105			END DO
106			CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
107			DO ij = 1, ip1jm
108			zv(ij) = ran1(idum) - .5
109			END DO
110			CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
111	guez	3
112	guez	26	CALL minmax(iip1*jjp1, zu, umin, ullm)
113			CALL minmax(iip1*jjm, zv, vmin, vllm)
114	guez	3
115	guez	26	ullm = abs(ullm)
116			vllm = abs(vllm)
117	guez	3
118	guez	26	DO l = 1, 50
119			IF (ii==1) THEN
120			IF (lstardis) THEN
121			CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
122			ELSE
123			CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
124			END IF
125			ELSE
126			IF (lstardis) THEN
127			CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
128			ELSE
129			CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
130			END IF
131			END IF
132	guez	3
133	guez	26	CALL minmax(iip1*jjp1, zu, umin, ullm)
134			CALL minmax(iip1*jjm, zv, vmin, vllm)
135	guez	3
136	guez	26	ullm = abs(ullm)
137			vllm = abs(vllm)
138	guez	3
139	guez	26	zllm = max(ullm, vllm)
140			z1llm = 1./zllm
141			DO ij = 1, ip1jmp1
142			zu(ij) = zu(ij)*z1llm
143			END DO
144			DO ij = 1, ip1jm
145			zv(ij) = zv(ij)*z1llm
146			END DO
147			end DO
148	guez	3
149	guez	26	IF (ii==1) THEN
150			IF (lstardis) THEN
151			cdivu = 1./zllm
152			ELSE
153			cdivu = zllm**(-1./nitergdiv)
154			END IF
155			ELSE
156			IF (lstardis) THEN
157			crot = 1./zllm
158			ELSE
159			crot = zllm**(-1./nitergrot)
160			END IF
161			END IF
162			END DO
163	guez	3
164	guez	26	PRINT *, 'cdivu = ', cdivu
165			PRINT *, 'crot = ', crot
166			PRINT *, 'cdivh = ', cdivh
167	guez	3
168	guez	26	! Variation verticale du coefficient de dissipation :
169	guez	27	zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
170			! (between 1 and 2)
171	guez	3
172	guez	26	tetaudiv = zvert / tetagdiv
173			tetaurot = zvert / tetagrot
174			tetah = zvert / tetatemp
175			call new_unit(unit)
176			open(unit, file="inidissip.csv", status="replace", action="write")
177			write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
178			do l = 1, llm
179			write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
180			end do
181			close(unit)
182			print *, 'Created file "inidissip.csv".'
183	guez	3
184	guez	27	max_zvert = maxval(zvert)
185			tetamin = min(1E6, tetagdiv / max_zvert, tetagrot / max_zvert, &
186			tetatemp / max_zvert)
187	guez	26	PRINT *, 'tetamin = ', tetamin
188	guez	27	idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
189	guez	26	PRINT *, 'idissip = ', idissip
190			dtdiss = idissip * dtvr
191			PRINT *, 'dtdiss = ', dtdiss
192
193			END SUBROUTINE inidissip
194
195			end module inidissip_m