1 | ! |
---|
2 | ! $Id: inidissip.F90 1707 2013-01-11 09:19:19Z dcugnet $ |
---|
3 | ! |
---|
4 | SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , & |
---|
5 | tetagdiv,tetagrot,tetatemp, vert_prof_dissip) |
---|
6 | !======================================================================= |
---|
7 | ! initialisation de la dissipation horizontale |
---|
8 | !======================================================================= |
---|
9 | !----------------------------------------------------------------------- |
---|
10 | ! declarations: |
---|
11 | ! ------------- |
---|
12 | |
---|
13 | USE control_mod, only : dissip_period,iperiod |
---|
14 | |
---|
15 | IMPLICIT NONE |
---|
16 | include "dimensions.h" |
---|
17 | include "paramet.h" |
---|
18 | include "comdissipn.h" |
---|
19 | include "comconst.h" |
---|
20 | include "comvert.h" |
---|
21 | include "logic.h" |
---|
22 | include "iniprint.h" |
---|
23 | |
---|
24 | LOGICAL,INTENT(in) :: lstardis |
---|
25 | INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh |
---|
26 | REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp |
---|
27 | |
---|
28 | integer, INTENT(in):: vert_prof_dissip |
---|
29 | ! Vertical profile of horizontal dissipation |
---|
30 | ! Allowed values: |
---|
31 | ! 0: rational fraction, function of pressure |
---|
32 | ! 1: tanh of altitude |
---|
33 | |
---|
34 | ! Local variables: |
---|
35 | REAL fact,zvert(llm),zz |
---|
36 | REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1) |
---|
37 | real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm) |
---|
38 | REAL ullm,vllm,umin,vmin,zhmin,zhmax |
---|
39 | REAL zllm |
---|
40 | |
---|
41 | INTEGER l,ij,idum,ii |
---|
42 | REAL tetamin |
---|
43 | REAL pseudoz |
---|
44 | character (len=80) :: abort_message |
---|
45 | |
---|
46 | REAL ran1 |
---|
47 | |
---|
48 | |
---|
49 | !----------------------------------------------------------------------- |
---|
50 | ! |
---|
51 | ! calcul des valeurs propres des operateurs par methode iterrative: |
---|
52 | ! ----------------------------------------------------------------- |
---|
53 | |
---|
54 | crot = -1. |
---|
55 | cdivu = -1. |
---|
56 | cdivh = -1. |
---|
57 | |
---|
58 | ! calcul de la valeur propre de divgrad: |
---|
59 | ! -------------------------------------- |
---|
60 | idum = 0 |
---|
61 | DO l = 1, llm |
---|
62 | DO ij = 1, ip1jmp1 |
---|
63 | deltap(ij,l) = 1. |
---|
64 | ENDDO |
---|
65 | ENDDO |
---|
66 | |
---|
67 | idum = -1 |
---|
68 | zh(1) = RAN1(idum)-.5 |
---|
69 | idum = 0 |
---|
70 | DO ij = 2, ip1jmp1 |
---|
71 | zh(ij) = RAN1(idum) -.5 |
---|
72 | ENDDO |
---|
73 | |
---|
74 | CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) |
---|
75 | |
---|
76 | CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) |
---|
77 | |
---|
78 | IF ( zhmin .GE. zhmax ) THEN |
---|
79 | write(lunout,*)' Inidissip zh min max ',zhmin,zhmax |
---|
80 | abort_message='probleme generateur alleatoire dans inidissip' |
---|
81 | call abort_gcm('inidissip',abort_message,1) |
---|
82 | ENDIF |
---|
83 | |
---|
84 | zllm = ABS( zhmax ) |
---|
85 | DO l = 1,50 |
---|
86 | IF(lstardis) THEN |
---|
87 | CALL divgrad2(1,zh,deltap,niterh,divgra) |
---|
88 | ELSE |
---|
89 | CALL divgrad (1,zh,niterh,divgra) |
---|
90 | ENDIF |
---|
91 | |
---|
92 | zllm = ABS(maxval(divgra)) |
---|
93 | zh = divgra / zllm |
---|
94 | ENDDO |
---|
95 | |
---|
96 | IF(lstardis) THEN |
---|
97 | cdivh = 1./ zllm |
---|
98 | ELSE |
---|
99 | cdivh = zllm ** ( -1./niterh ) |
---|
100 | ENDIF |
---|
101 | |
---|
102 | ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) |
---|
103 | ! ----------------------------------------------------------------- |
---|
104 | write(lunout,*)'inidissip: calcul des valeurs propres' |
---|
105 | |
---|
106 | DO ii = 1, 2 |
---|
107 | ! |
---|
108 | DO ij = 1, ip1jmp1 |
---|
109 | zu(ij) = RAN1(idum) -.5 |
---|
110 | ENDDO |
---|
111 | CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) |
---|
112 | DO ij = 1, ip1jm |
---|
113 | zv(ij) = RAN1(idum) -.5 |
---|
114 | ENDDO |
---|
115 | CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) |
---|
116 | |
---|
117 | CALL minmax(iip1*jjp1,zu,umin,ullm ) |
---|
118 | CALL minmax(iip1*jjm, zv,vmin,vllm ) |
---|
119 | |
---|
120 | ullm = ABS ( ullm ) |
---|
121 | vllm = ABS ( vllm ) |
---|
122 | |
---|
123 | DO l = 1, 50 |
---|
124 | IF(ii.EQ.1) THEN |
---|
125 | !cccc CALL covcont( 1,zu,zv,zu,zv ) |
---|
126 | IF(lstardis) THEN |
---|
127 | CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy ) |
---|
128 | ELSE |
---|
129 | CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy ) |
---|
130 | ENDIF |
---|
131 | ELSE |
---|
132 | IF(lstardis) THEN |
---|
133 | CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy ) |
---|
134 | ELSE |
---|
135 | CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy ) |
---|
136 | ENDIF |
---|
137 | ENDIF |
---|
138 | |
---|
139 | zllm = max(abs(maxval(gx)), abs(maxval(gy))) |
---|
140 | zu = gx / zllm |
---|
141 | zv = gy / zllm |
---|
142 | end DO |
---|
143 | |
---|
144 | IF ( ii.EQ.1 ) THEN |
---|
145 | IF(lstardis) THEN |
---|
146 | cdivu = 1./zllm |
---|
147 | ELSE |
---|
148 | cdivu = zllm **( -1./nitergdiv ) |
---|
149 | ENDIF |
---|
150 | ELSE |
---|
151 | IF(lstardis) THEN |
---|
152 | crot = 1./ zllm |
---|
153 | ELSE |
---|
154 | crot = zllm **( -1./nitergrot ) |
---|
155 | ENDIF |
---|
156 | ENDIF |
---|
157 | |
---|
158 | end DO |
---|
159 | |
---|
160 | ! petit test pour les operateurs non star: |
---|
161 | ! ---------------------------------------- |
---|
162 | |
---|
163 | ! IF(.NOT.lstardis) THEN |
---|
164 | fact = rad*24./REAL(jjm) |
---|
165 | fact = fact*fact |
---|
166 | write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu |
---|
167 | write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot |
---|
168 | write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh |
---|
169 | ! ENDIF |
---|
170 | |
---|
171 | !----------------------------------------------------------------------- |
---|
172 | ! variation verticale du coefficient de dissipation: |
---|
173 | ! -------------------------------------------------- |
---|
174 | |
---|
175 | if (vert_prof_dissip == 1) then |
---|
176 | do l=1,llm |
---|
177 | pseudoz=8.*log(preff/presnivs(l)) |
---|
178 | zvert(l)=1+ & |
---|
179 | (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. & |
---|
180 | *(dissip_factz-1.) |
---|
181 | enddo |
---|
182 | else |
---|
183 | DO l=1,llm |
---|
184 | zvert(l)=1. |
---|
185 | ENDDO |
---|
186 | fact=2. |
---|
187 | DO l = 1, llm |
---|
188 | zz = 1. - preff/presnivs(l) |
---|
189 | zvert(l)= fact -( fact-1.)/( 1.+zz*zz ) |
---|
190 | ENDDO |
---|
191 | endif |
---|
192 | |
---|
193 | |
---|
194 | write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' |
---|
195 | |
---|
196 | tetamin = 1.e+6 |
---|
197 | |
---|
198 | DO l=1,llm |
---|
199 | tetaudiv(l) = zvert(l)/tetagdiv |
---|
200 | tetaurot(l) = zvert(l)/tetagrot |
---|
201 | tetah(l) = zvert(l)/tetatemp |
---|
202 | |
---|
203 | IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) |
---|
204 | IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) |
---|
205 | IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) |
---|
206 | ENDDO |
---|
207 | |
---|
208 | ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def |
---|
209 | IF (dissip_period == 0) THEN |
---|
210 | dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod |
---|
211 | write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period |
---|
212 | dissip_period = MAX(iperiod,dissip_period) |
---|
213 | END IF |
---|
214 | |
---|
215 | dtdiss = dissip_period * dtvr |
---|
216 | write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr |
---|
217 | |
---|
218 | DO l = 1,llm |
---|
219 | write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), & |
---|
220 | dtdiss*tetah(l) |
---|
221 | ENDDO |
---|
222 | |
---|
223 | END SUBROUTINE inidissip |
---|