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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 54 - (hide 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 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 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     USE paramet_m, ONLY : jjp1
30 guez 48 use jumble, only: new_unit
31 guez 27 use filtreg_m, only: filtreg
32 guez 54 use gradiv2_m, only: gradiv2
33 guez 3
34 guez 26 ! Variables local to the procedure:
35 guez 27 REAL zvert(llm), max_zvert
36 guez 54 REAL, dimension(iim + 1, jjm + 1):: zh, zu
37     real zv(iim + 1, jjm), deltap(iim + 1, jjm + 1, llm)
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 54 CALL divgrad2(1, zh, deltap, niterh, zh, -1.)
57 guez 26 ELSE
58 guez 54 CALL divgrad(1, zh, niterh, zh, -1.)
59 guez 26 END IF
60 guez 3
61 guez 40 zllm = abs(maxval(zh))
62     zh = zh / 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     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 guez 3
87 guez 54 zllm = max(abs(maxval(zu)), abs(maxval(zv)))
88     zu = zu / zllm
89     zv = zv / zllm
90     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     CALL nxgraro2(1, zu, zv, nitergrot, zu, zv, -1.)
110 guez 26 ELSE
111 guez 54 CALL nxgrarot(1, zu, zv, nitergrot, zu, zv, -1.)
112 guez 26 END IF
113 guez 3
114 guez 54 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 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     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 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     PRINT *, 'dtdiss = ', dtdiss
151    
152     END SUBROUTINE inidissip
153    
154     end module inidissip_m
  ViewVC Help
Powered by ViewVC 1.1.21