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 *, '  Inidissip  zh min max  ', zhmin, zhmax
       STOP 'probleme generateur alleatoire dans inidissip'
    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	3
70	guez	26	IF (zhmin>=zhmax) THEN
71			PRINT *, ' Inidissip zh min max ', zhmin, zhmax
72			STOP 'probleme generateur alleatoire dans inidissip'
73			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	guez	3
104	guez	26	DO ij = 1, ip1jmp1
105			zu(ij) = ran1(idum) - .5
106			END DO
107			CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
108			DO ij = 1, ip1jm
109			zv(ij) = ran1(idum) - .5
110			END DO
111			CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
112	guez	3
113	guez	26	CALL minmax(iip1*jjp1, zu, umin, ullm)
114			CALL minmax(iip1*jjm, zv, vmin, vllm)
115	guez	3
116	guez	26	ullm = abs(ullm)
117			vllm = abs(vllm)
118	guez	3
119	guez	26	DO l = 1, 50
120			IF (ii==1) THEN
121			IF (lstardis) THEN
122			CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
123			ELSE
124			CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
125			END IF
126			ELSE
127			IF (lstardis) THEN
128			CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
129			ELSE
130			CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
131			END IF
132			END IF
133	guez	3
134	guez	26	CALL minmax(iip1*jjp1, zu, umin, ullm)
135			CALL minmax(iip1*jjm, zv, vmin, vllm)
136	guez	3
137	guez	26	ullm = abs(ullm)
138			vllm = abs(vllm)
139	guez	3
140	guez	26	zllm = max(ullm, vllm)
141			z1llm = 1./zllm
142			DO ij = 1, ip1jmp1
143			zu(ij) = zu(ij)*z1llm
144			END DO
145			DO ij = 1, ip1jm
146			zv(ij) = zv(ij)*z1llm
147			END DO
148			end DO
149	guez	3
150	guez	26	IF (ii==1) THEN
151			IF (lstardis) THEN
152			cdivu = 1./zllm
153			ELSE
154			cdivu = zllm**(-1./nitergdiv)
155			END IF
156			ELSE
157			IF (lstardis) THEN
158			crot = 1./zllm
159			ELSE
160			crot = zllm**(-1./nitergrot)
161			END IF
162			END IF
163	guez	3
164	guez	26	END DO
165	guez	3
166	guez	26	PRINT *, 'cdivu = ', cdivu
167			PRINT *, 'crot = ', crot
168			PRINT *, 'cdivh = ', cdivh
169	guez	3
170	guez	26	! Variation verticale du coefficient de dissipation :
171	guez	27	zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
172			! (between 1 and 2)
173	guez	3
174	guez	26	tetaudiv = zvert / tetagdiv
175			tetaurot = zvert / tetagrot
176			tetah = zvert / tetatemp
177			call new_unit(unit)
178			open(unit, file="inidissip.csv", status="replace", action="write")
179			write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
180			do l = 1, llm
181			write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
182			end do
183			close(unit)
184			print *, 'Created file "inidissip.csv".'
185	guez	3
186	guez	27	max_zvert = maxval(zvert)
187			tetamin = min(1E6, tetagdiv / max_zvert, tetagrot / max_zvert, &
188			tetatemp / max_zvert)
189	guez	26	PRINT *, 'tetamin = ', tetamin
190	guez	27	idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
191	guez	26	PRINT *, 'idissip = ', idissip
192			dtdiss = idissip * dtvr
193			PRINT *, 'dtdiss = ', dtdiss
194
195			END SUBROUTINE inidissip
196
197			end module inidissip_m