dyn3d/Dissipation/inidissip.f

module inidissip_m

  use dimens_m, only: llm

  IMPLICIT NONE

  private llm

  REAL dtdiss ! in s
  integer idissip ! période de la dissipation (en pas de temps)
  real tetaudiv(llm), tetaurot(llm), tetah(llm) ! in s
  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. Calcul des valeurs
    ! propres des opérateurs par méthode itérative.

    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 : iim, jjm, llm
    use filtreg_m, only: filtreg
    use gradiv2_m, only: gradiv2
    use jumble, only: new_unit
    USE paramet_m, ONLY : jjp1

    ! Variables local to the procedure:
    REAL zvert(llm), max_zvert ! no dimension
    REAL, dimension(iim + 1, jjm + 1, 1):: zh, zu, gx, divgra, deltap
    real zv(iim + 1, jjm, 1), gy(iim + 1, jjm, 1)
    REAL zllm
    INTEGER l, seed_size, ii, unit
    REAL tetamin ! in s

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

    PRINT *, 'Call sequence information: inidissip'
    call random_seed(size=seed_size)
    call random_seed(put=(/(0, ii = 1, seed_size)/))

    PRINT *, 'Calcul des valeurs propres de divgrad'
    deltap = 1.
    call random_number(zh)
    zh = zh - 0.5
    CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)

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

       zllm = abs(maxval(divgra))
       zh = divgra / zllm
    END DO

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

    PRINT *, 'Calcul des valeurs propres de gradiv'
    call random_number(zu)
    zu = zu - 0.5
    CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
    call random_number(zv)
    zv = zv - 0.5
    CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)

    DO l = 1, 50
       IF (lstardis) THEN
          CALL gradiv2(zu, zv, nitergdiv, gx, gy, -1.)
       ELSE
          CALL gradiv(1, zu, zv, nitergdiv, gx, gy, -1.)
       END IF

       zllm = max(abs(maxval(gx)), abs(maxval(gy)))
       zu = gx / zllm
       zv = gy / zllm
    end DO

    IF (lstardis) THEN
       cdivu = 1. / zllm
    ELSE
       cdivu = zllm**(- 1. / nitergdiv)
    END IF
    PRINT *, 'cdivu = ', cdivu

    PRINT *, 'Calcul des valeurs propres de nxgrarot'
    call random_number(zu)
    zu = zu - 0.5
    CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
    call random_number(zv)
    zv = zv - 0.5
    CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)

    DO l = 1, 50
       IF (lstardis) THEN
          CALL nxgraro2(1, zu, zv, nitergrot, gx, gy, -1.)
       ELSE
          CALL nxgrarot(1, zu, zv, nitergrot, gx, gy, -1.)
       END IF

       zllm = max(abs(maxval(gx)), abs(maxval(gy)))
       zu = gx / zllm
       zv = gy / zllm
    end DO

    IF (lstardis) THEN
       crot = 1. / zllm
    ELSE
       crot = zllm**(- 1. / nitergrot)
    END IF
    PRINT *, 'crot = ', crot

    ! 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 (s)" "tetaurot (s)" "tetah (s)"' ! 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, "s"

  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	57	REAL dtdiss ! in s
10	guez	47	integer idissip ! période de la dissipation (en pas de temps)
11	guez	57	real tetaudiv(llm), tetaurot(llm), tetah(llm) ! in s
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	3
20	guez	54	! Initialisation de la dissipation horizontale. Calcul des valeurs
21			! propres des opérateurs par méthode itérative.
22
23	guez	26	USE comconst, ONLY : dtvr
24	guez	27	use comdissnew, only: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
25			tetagrot, tetatemp
26	guez	26	USE comvert, ONLY : preff, presnivs
27			USE conf_gcm_m, ONLY : iperiod
28	guez	54	USE dimens_m, ONLY : iim, jjm, llm
29	guez	27	use filtreg_m, only: filtreg
30	guez	54	use gradiv2_m, only: gradiv2
31	guez	57	use jumble, only: new_unit
32			USE paramet_m, ONLY : jjp1
33	guez	3
34	guez	26	! Variables local to the procedure:
35	guez	57	REAL zvert(llm), max_zvert ! no dimension
36			REAL, dimension(iim + 1, jjm + 1, 1):: zh, zu, gx, divgra, deltap
37			real zv(iim + 1, jjm, 1), gy(iim + 1, jjm, 1)
38	guez	40	REAL zllm
39	guez	54	INTEGER l, seed_size, ii, unit
40	guez	27	REAL tetamin ! in s
41	guez	3
42	guez	26	!-----------------------------------------------------------------------
43	guez	3
44	guez	26	PRINT *, 'Call sequence information: inidissip'
45	guez	54	call random_seed(size=seed_size)
46			call random_seed(put=(/(0, ii = 1, seed_size)/))
47	guez	3
48	guez	54	PRINT *, 'Calcul des valeurs propres de divgrad'
49	guez	40	deltap = 1.
50	guez	54	call random_number(zh)
51			zh = zh - 0.5
52	guez	26	CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
53	guez	3
54	guez	26	DO l = 1, 50
55			IF (lstardis) THEN
56	guez	57	CALL divgrad2(1, zh, deltap, niterh, divgra, -1.)
57	guez	26	ELSE
58	guez	57	CALL divgrad(1, zh, niterh, divgra, -1.)
59	guez	26	END IF
60	guez	3
61	guez	57	zllm = abs(maxval(divgra))
62			zh = divgra / zllm
63	guez	26	END DO
64	guez	3
65	guez	26	IF (lstardis) THEN
66	guez	40	cdivh = 1. / zllm
67	guez	26	ELSE
68	guez	40	cdivh = zllm**(- 1. / niterh)
69	guez	26	END IF
70	guez	54	PRINT *, 'cdivh = ', cdivh
71	guez	3
72	guez	54	PRINT *, 'Calcul des valeurs propres de gradiv'
73			call random_number(zu)
74			zu = zu - 0.5
75			CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
76			call random_number(zv)
77			zv = zv - 0.5
78			CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
79	guez	3
80	guez	54	DO l = 1, 50
81			IF (lstardis) THEN
82	guez	57	CALL gradiv2(zu, zv, nitergdiv, gx, gy, -1.)
83	guez	54	ELSE
84	guez	57	CALL gradiv(1, zu, zv, nitergdiv, gx, gy, -1.)
85	guez	54	END IF
86	guez	3
87	guez	57	zllm = max(abs(maxval(gx)), abs(maxval(gy)))
88			zu = gx / zllm
89			zv = gy / zllm
90	guez	54	end DO
91	guez	3
92	guez	54	IF (lstardis) THEN
93			cdivu = 1. / zllm
94			ELSE
95			cdivu = zllm**(- 1. / nitergdiv)
96			END IF
97			PRINT *, 'cdivu = ', cdivu
98	guez	3
99	guez	54	PRINT *, 'Calcul des valeurs propres de nxgrarot'
100			call random_number(zu)
101			zu = zu - 0.5
102			CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1)
103			call random_number(zv)
104			zv = zv - 0.5
105			CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1)
106	guez	3
107	guez	54	DO l = 1, 50
108			IF (lstardis) THEN
109	guez	57	CALL nxgraro2(1, zu, zv, nitergrot, gx, gy, -1.)
110	guez	26	ELSE
111	guez	57	CALL nxgrarot(1, zu, zv, nitergrot, gx, gy, -1.)
112	guez	26	END IF
113	guez	3
114	guez	57	zllm = max(abs(maxval(gx)), abs(maxval(gy)))
115			zu = gx / zllm
116			zv = gy / zllm
117	guez	54	end DO
118
119			IF (lstardis) THEN
120			crot = 1. / zllm
121			ELSE
122	guez	57	crot = zllm**(- 1. / nitergrot)
123	guez	54	END IF
124	guez	26	PRINT *, 'crot = ', crot
125	guez	3
126	guez	26	! Variation verticale du coefficient de dissipation :
127	guez	27	zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
128			! (between 1 and 2)
129	guez	3
130	guez	26	tetaudiv = zvert / tetagdiv
131			tetaurot = zvert / tetagrot
132			tetah = zvert / tetatemp
133	guez	54
134	guez	26	call new_unit(unit)
135			open(unit, file="inidissip.csv", status="replace", action="write")
136	guez	57	write(unit, fmt=*) '"tetaudiv (s)" "tetaurot (s)" "tetah (s)"' ! title line
137	guez	26	do l = 1, llm
138			write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l)
139			end do
140			close(unit)
141			print *, 'Created file "inidissip.csv".'
142	guez	3
143	guez	27	max_zvert = maxval(zvert)
144	guez	54	tetamin = min(1e6, tetagdiv / max_zvert, tetagrot / max_zvert, &
145	guez	27	tetatemp / max_zvert)
146	guez	26	PRINT *, 'tetamin = ', tetamin
147	guez	27	idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
148	guez	26	PRINT *, 'idissip = ', idissip
149			dtdiss = idissip * dtvr
150	guez	57	PRINT *, 'dtdiss = ', dtdiss, "s"
151	guez	26
152			END SUBROUTINE inidissip
153
154			end module inidissip_m