dyn3d/Dissipation/inidissip.f90

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-1
  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: nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, &
         tetatemp
    USE disvert_m, ONLY: preff, presnivs
    USE conf_gcm_m, ONLY: iperiod
    USE dimens_m, ONLY: iim, jjm
    use divgrad2_m, only: divgrad2
    use filtreg_m, only: filtreg
    use gradiv2_m, only: gradiv2
    use jumble, only: new_unit
    use nxgraro2_m, only: nxgraro2
    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=(/(1, 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.)

    DO l = 1, 50
       CALL divgrad2(1, zh, deltap, niterh, divgra, -1.)
       zllm = abs(maxval(divgra))
       zh = divgra / zllm
    END DO

    cdivh = 1. / zllm
    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.)
    call random_number(zv)
    zv = zv - 0.5
    CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.)

    DO l = 1, 50
       CALL gradiv2(zu, zv, nitergdiv, gx, gy, -1.)
       zllm = max(abs(maxval(gx)), abs(maxval(gy)))
       zu = gx / zllm
       zv = gy / zllm
    end DO

    cdivu = 1. / zllm
    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.)
    call random_number(zv)
    zv = zv - 0.5
    CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.)

    DO l = 1, 50
       CALL nxgraro2(zu, zv, nitergrot, gx, gy, -1.)
       zllm = max(abs(maxval(gx)), abs(maxval(gy)))
       zu = gx / zllm
       zv = gy / zllm
    end DO

    crot = 1. / zllm
    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

    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"

    call new_unit(unit)
    open(unit, file="inidissip.csv", status="replace", action="write")

    ! Title line:
    write(unit, fmt=*) '"presnivs (hPa)" "dtdiss * tetaudiv" ' &
         // '"dtdiss * tetaurot" "dtdiss * tetah"'

    do l = 1, llm
       write(unit, fmt=*) presnivs(l) / 100., dtdiss * tetaudiv(l), &
            dtdiss * tetaurot(l), dtdiss * tetah(l)
    end do
    close(unit)
    print *, 'Created file "inidissip.csv".'

  END SUBROUTINE inidissip

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