dyn3d/Dissipation/inidissip.f

module inidissip_m

  use dimens_m, only: llm

  IMPLICIT NONE

  private llm

  REAL dtdiss
  integer idissip ! période de la dissipation (en pas de temps)
  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 jumble, 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 zhmin, zhmax
    REAL zllm
    INTEGER l, ij, idum, ii, unit
    REAL tetamin ! in s
    REAL ran1

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

    PRINT *, 'Call sequence information: inidissip'

    ! Calcul des valeurs propres des opérateurs par méthode itérative :

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

    ! Calcul de la valeur propre de divgrad :

    deltap = 1.
    idum = -1
    zh(1) = ran1(idum) - 0.5
    idum = 0
    DO ij = 2, ip1jmp1
       zh(ij) = ran1(idum) - 0.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

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

       zllm = abs(maxval(zh))
       zh = zh / zllm
    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) - 0.5
       END DO
       CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
       DO ij = 1, ip1jm
          zv(ij) = ran1(idum) - 0.5
       END DO
       CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)

       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

          zllm = max(abs(maxval(zu)), abs(maxval(zv)))
          zu = zu / zllm
          zv = zv / zllm
       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	guez	47	integer idissip ! période de la dissipation (en pas de temps)
11	guez	40	real tetaudiv(llm), tetaurot(llm), tetah(llm)
12	guez	26	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	guez	40	! 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	guez	48	use jumble, 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	guez	40	REAL zhmin, zhmax
35			REAL zllm
36	guez	26	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	40	! Calcul des valeurs propres des opérateurs par méthode itérative :
45	guez	3
46	guez	26	crot = -1.
47			cdivu = -1.
48			cdivh = -1.
49	guez	3
50	guez	40	! Calcul de la valeur propre de divgrad :
51	guez	3
52	guez	40	deltap = 1.
53	guez	26	idum = -1
54	guez	40	zh(1) = ran1(idum) - 0.5
55	guez	26	idum = 0
56			DO ij = 2, ip1jmp1
57	guez	40	zh(ij) = ran1(idum) - 0.5
58	guez	26	END DO
59	guez	3
60	guez	26	CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
61	guez	3
62	guez	26	CALL minmax(iip1*jjp1, zh, zhmin, zhmax)
63	guez	39	IF (zhmin >= zhmax) THEN
64			PRINT *, 'zhmin zhmax', zhmin, zhmax
65			print *, 'Problème générateur aléatoire dans inidissip'
66			STOP 1
67	guez	26	END IF
68	guez	3
69	guez	26	DO l = 1, 50
70			IF (lstardis) THEN
71			CALL divgrad2(1, zh, deltap, niterh, zh)
72			ELSE
73			CALL divgrad(1, zh, niterh, zh)
74			END IF
75	guez	3
76	guez	40	zllm = abs(maxval(zh))
77			zh = zh / zllm
78	guez	26	END DO
79	guez	3
80	guez	26	IF (lstardis) THEN
81	guez	40	cdivh = 1. / zllm
82	guez	26	ELSE
83	guez	40	cdivh = zllm**(- 1. / niterh)
84	guez	26	END IF
85	guez	3
86	guez	40	! Calcul des valeurs propres de gradiv (ii = 1) et nxgrarot (ii = 2)
87	guez	3
88	guez	40	PRINT *, 'Calcul des valeurs propres'
89	guez	3
90	guez	40	DO ii = 1, 2
91	guez	26	DO ij = 1, ip1jmp1
92	guez	40	zu(ij) = ran1(idum) - 0.5
93	guez	26	END DO
94			CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
95			DO ij = 1, ip1jm
96	guez	40	zv(ij) = ran1(idum) - 0.5
97	guez	26	END DO
98			CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
99	guez	3
100	guez	40	DO l = 1, 50
101	guez	26	IF (ii==1) THEN
102			IF (lstardis) THEN
103			CALL gradiv2(1, zu, zv, nitergdiv, zu, zv)
104			ELSE
105			CALL gradiv(1, zu, zv, nitergdiv, zu, zv)
106			END IF
107			ELSE
108			IF (lstardis) THEN
109			CALL nxgraro2(1, zu, zv, nitergrot, zu, zv)
110			ELSE
111			CALL nxgrarot(1, zu, zv, nitergrot, zu, zv)
112			END IF
113			END IF
114	guez	3
115	guez	40	zllm = max(abs(maxval(zu)), abs(maxval(zv)))
116			zu = zu / zllm
117			zv = zv / zllm
118	guez	26	end DO
119	guez	3
120	guez	26	IF (ii==1) THEN
121			IF (lstardis) THEN
122	guez	40	cdivu = 1. / zllm
123	guez	26	ELSE
124	guez	40	cdivu = zllm**(- 1. / nitergdiv)
125	guez	26	END IF
126			ELSE
127			IF (lstardis) THEN
128			crot = 1./zllm
129			ELSE
130	guez	40	crot = zllm**(-1. / nitergrot)
131	guez	26	END IF
132			END IF
133			END DO
134	guez	3
135	guez	26	PRINT *, 'cdivu = ', cdivu
136			PRINT *, 'crot = ', crot
137			PRINT *, 'cdivh = ', cdivh
138	guez	3
139	guez	26	! Variation verticale du coefficient de dissipation :
140	guez	27	zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
141			! (between 1 and 2)
142	guez	3
143	guez	26	tetaudiv = zvert / tetagdiv
144			tetaurot = zvert / tetagrot
145			tetah = zvert / tetatemp
146			call new_unit(unit)
147			open(unit, file="inidissip.csv", status="replace", action="write")
148			write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
149			do l = 1, llm
150			write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
151			end do
152			close(unit)
153			print *, 'Created file "inidissip.csv".'
154	guez	3
155	guez	27	max_zvert = maxval(zvert)
156			tetamin = min(1E6, tetagdiv / max_zvert, tetagrot / max_zvert, &
157			tetatemp / max_zvert)
158	guez	26	PRINT *, 'tetamin = ', tetamin
159	guez	27	idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
160	guez	26	PRINT *, 'idissip = ', idissip
161			dtdiss = idissip * dtvr
162			PRINT *, 'dtdiss = ', dtdiss
163
164			END SUBROUTINE inidissip
165
166			end module inidissip_m