/[lmdze]/trunk/dyn3d/Dissipation/inidissip.f
ViewVC logotype

Contents of /trunk/dyn3d/Dissipation/inidissip.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 54 - (show annotations)
Tue Dec 6 15:07:04 2011 UTC (12 years, 5 months ago) by guez
Original Path: trunk/libf/dyn3d/inidissip.f90
File size: 4110 byte(s) 
Removed Numerical Recipes procedure "ran1". Replaced calls to "ran1"
in "inidissip" by calls to intrinsic procedures.

Split file "interface_surf.f90" into a file with a module containing
only variables, "interface_surf", and single-procedure files. Gathered
files into directory "Interface_surf".

Added argument "cdivu" to "gradiv" and "gradiv2", "cdivh" to
"divgrad2" and "divgrad", and "crot" to "nxgraro2" and
"nxgrarot". "dissip" now uses variables "cdivu", "cdivh" and "crot"
from module "inidissip_m", so it can pass them to "gradiv2",
etc. Thanks to this modification, we avoid a circular dependency
betwwen "inidissip.f90" and "gradiv2.f90", etc. The value -1. used by
"gradiv2", for instance, during computation of eigenvalues is not the
value "cdivu" computed by "inidissip".

Extracted procedure "start_inter_3d" from module "startdyn", to its
own module.

In "inidissip", unrolled loop on "ii". I find it clearer now.

Moved variables "matriceun", "matriceus", "matricevn", "matricevs",
"matrinvn" and "matrinvs" from module "parafilt" to module
"inifilr_m". Moved variables "jfiltnu", "jfiltnv", "jfiltsu",
"jfiltsv" from module "coefils" to module "inifilr_m".
1 module inidissip_m
2
3 use dimens_m, only: llm
4
5 IMPLICIT NONE
6
7 private llm
8
9 REAL dtdiss
10 integer idissip ! période de la dissipation (en pas de temps)
11 real tetaudiv(llm), tetaurot(llm), tetah(llm)
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: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
25 tetagrot, tetatemp
26 USE comvert, ONLY : preff, presnivs
27 USE conf_gcm_m, ONLY : iperiod
28 USE dimens_m, ONLY : iim, jjm, llm
29 USE paramet_m, ONLY : jjp1
30 use jumble, only: new_unit
31 use filtreg_m, only: filtreg
32 use gradiv2_m, only: gradiv2
33
34 ! Variables local to the procedure:
35 REAL zvert(llm), max_zvert
36 REAL, dimension(iim + 1, jjm + 1):: zh, zu
37 real zv(iim + 1, jjm), deltap(iim + 1, jjm + 1, llm)
38 REAL zllm
39 INTEGER l, seed_size, ii, unit
40 REAL tetamin ! in s
41
42 !-----------------------------------------------------------------------
43
44 PRINT *, 'Call sequence information: inidissip'
45 call random_seed(size=seed_size)
46 call random_seed(put=(/(0, ii = 1, seed_size)/))
47
48 PRINT *, 'Calcul des valeurs propres de divgrad'
49 deltap = 1.
50 call random_number(zh)
51 zh = zh - 0.5
52 CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1)
53
54 DO l = 1, 50
55 IF (lstardis) THEN
56 CALL divgrad2(1, zh, deltap, niterh, zh, -1.)
57 ELSE
58 CALL divgrad(1, zh, niterh, zh, -1.)
59 END IF
60
61 zllm = abs(maxval(zh))
62 zh = zh / zllm
63 END DO
64
65 IF (lstardis) THEN
66 cdivh = 1. / zllm
67 ELSE
68 cdivh = zllm**(- 1. / niterh)
69 END IF
70 PRINT *, 'cdivh = ', cdivh
71
72 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
80 DO l = 1, 50
81 IF (lstardis) THEN
82 CALL gradiv2(1, zu, zv, nitergdiv, zu, zv, -1.)
83 ELSE
84 CALL gradiv(1, zu, zv, nitergdiv, zu, zv, -1.)
85 END IF
86
87 zllm = max(abs(maxval(zu)), abs(maxval(zv)))
88 zu = zu / zllm
89 zv = zv / zllm
90 end DO
91
92 IF (lstardis) THEN
93 cdivu = 1. / zllm
94 ELSE
95 cdivu = zllm**(- 1. / nitergdiv)
96 END IF
97 PRINT *, 'cdivu = ', cdivu
98
99 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
107 DO l = 1, 50
108 IF (lstardis) THEN
109 CALL nxgraro2(1, zu, zv, nitergrot, zu, zv, -1.)
110 ELSE
111 CALL nxgrarot(1, zu, zv, nitergrot, zu, zv, -1.)
112 END IF
113
114 zllm = max(abs(maxval(zu)), abs(maxval(zv)))
115 zu = zu / zllm
116 zv = zv / zllm
117 end DO
118
119 IF (lstardis) THEN
120 crot = 1. / zllm
121 ELSE
122 crot = zllm**(-1. / nitergrot)
123 END IF
124 PRINT *, 'crot = ', crot
125
126 ! Variation verticale du coefficient de dissipation :
127 zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2)
128 ! (between 1 and 2)
129
130 tetaudiv = zvert / tetagdiv
131 tetaurot = zvert / tetagrot
132 tetah = zvert / tetatemp
133
134 call new_unit(unit)
135 open(unit, file="inidissip.csv", status="replace", action="write")
136 write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line
137 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
143 max_zvert = maxval(zvert)
144 tetamin = min(1e6, tetagdiv / max_zvert, tetagrot / max_zvert, &
145 tetatemp / max_zvert)
146 PRINT *, 'tetamin = ', tetamin
147 idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod
148 PRINT *, 'idissip = ', idissip
149 dtdiss = idissip * dtvr
150 PRINT *, 'dtdiss = ', dtdiss
151
152 END SUBROUTINE inidissip
153
154 end module inidissip_m
  ViewVC Help
Powered by ViewVC 1.1.21