/[lmdze]/trunk/libf/dyn3d/guide.f90
ViewVC logotype

Contents of /trunk/libf/dyn3d/guide.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 44 - (show annotations)
Wed Apr 13 12:29:18 2011 UTC (13 years, 1 month ago) by guez
File size: 10115 byte(s)
Removed argument "pdteta" of "calfis", because it was not used.

Created module "conf_guide_m", containing procedure
"conf_guide". Moved module variables from "guide_m" to "conf_guide_m".

In module "getparam", removed "ini_getparam" and "fin_getparam" from
generic interface "getpar".

Created module variables in "tau2alpha_m" to replace common "comdxdy".

1 MODULE guide_m
2
3 ! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09
4 ! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06
5
6 IMPLICIT NONE
7
8 REAL aire_min, aire_max
9
10 CONTAINS
11
12 SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)
13
14 ! Author: F.Hourdin
15
16 USE comconst, ONLY : cpp, daysec, dtvr, kappa
17 USE comgeom, ONLY : aire, rlatu, rlonv
18 USE comvert, ONLY : ap, bp, preff, presnivs
19 USE conf_gcm_m, ONLY : day_step, iperiod
20 use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &
21 guide_p, ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &
22 tau_min_t, tau_max_t, tau_min_q, tau_max_q, tau_min_p, tau_max_p, &
23 online
24 USE dimens_m, ONLY : jjm, llm
25 USE exner_hyb_m, ONLY : exner_hyb
26 USE inigrads_m, ONLY : inigrads
27 use netcdf, only: nf90_nowrite, nf90_open, nf90_close, nf90_inq_dimid, &
28 nf90_inquire_dimension
29 use nr_util, only: pi
30 USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
31 USE q_sat_m, ONLY : q_sat
32 USE serre, ONLY : clat, clon
33 use tau2alpha_m, only: tau2alpha, dxdys
34
35 ! variables dynamiques
36 REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
37 REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle
38 REAL q(ip1jmp1, llm) ! temperature potentielle
39 REAL ps(ip1jmp1) ! pression au sol
40 REAL masse(ip1jmp1, llm) ! masse d'air
41
42 ! variables dynamiques pour les reanalyses.
43 REAL, save:: ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas
44 REAL, save:: tetarea1(ip1jmp1, llm) ! temp pot reales
45 REAL, save:: qrea1(ip1jmp1, llm) ! temp pot reales
46 REAL, save:: psrea1(ip1jmp1) ! ps
47 REAL, save:: ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas
48 REAL, save:: tetarea2(ip1jmp1, llm) ! temp pot reales
49 REAL, save:: qrea2(ip1jmp1, llm) ! temp pot reales
50 REAL, save:: masserea2(ip1jmp1, llm) ! masse
51 REAL, save:: psrea2(ip1jmp1) ! ps
52
53 REAL, save:: alpha_q(ip1jmp1)
54 REAL, save:: alpha_t(ip1jmp1), alpha_p(ip1jmp1)
55 REAL, save:: alpha_u(ip1jmp1), alpha_v(ip1jm)
56 REAL dday_step, toto, reste
57 real, save:: itau_test
58 INTEGER, save:: step_rea, count_no_rea
59
60 INTEGER ilon, ilat
61 REAL factt, ztau(ip1jmp1)
62
63 INTEGER, INTENT(IN):: itau
64 INTEGER ij, l
65 INTEGER ncidpl, varidpl, status
66 INTEGER rcod, rid
67 REAL ditau, tau, a
68 INTEGER, SAVE:: nlev
69
70 ! TEST SUR QSAT
71 REAL p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)
72 REAL pkf(ip1jmp1, llm)
73 REAL pres(ip1jmp1, llm)
74
75 REAL qsat(ip1jmp1, llm)
76 REAL unskap
77 REAL tnat(ip1jmp1, llm)
78
79 LOGICAL:: first = .TRUE.
80 CHARACTER(len=10) file
81 INTEGER:: igrads = 2
82 REAL:: dtgrads = 100.
83
84 !-----------------------------------------------------------------------
85
86 PRINT *, 'Call sequence information: guide'
87
88 ! calcul de l'humidite saturante
89
90 forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
91 CALL massdair(p, masse)
92 CALL exner_hyb(ps, p, pks, pk, pkf)
93 tnat(:, :) = pk(:, :)*teta(:, :)/cpp
94 unskap = 1./kappa
95 pres(:, :) = preff*(pk(:, :)/cpp)**unskap
96 qsat = q_sat(tnat, pres)
97
98 ! initialisations pour la lecture des reanalyses.
99 ! alpha determine la part des injections de donnees a chaque etape
100 ! alpha=1 signifie pas d'injection
101 ! alpha=0 signifie injection totale
102
103 IF (online==-1) THEN
104 RETURN
105 END IF
106
107 IF (first) THEN
108 CALL conf_guide
109 file = 'guide'
110 CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
111 180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
112 PRINT *, '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
113 IF (online==-1) RETURN
114
115 IF (online==1) THEN
116 ! Constantes de temps de rappel en jour
117 ! 0.1 c'est en gros 2h30.
118 ! 1e10 est une constante infinie donc en gros pas de guidage
119
120 ! coordonnees du centre du zoom
121 CALL coordij(clon, clat, ilon, ilat)
122 ! aire de la maille au centre du zoom
123 aire_min = aire(ilon+(ilat-1)*iip1)
124 ! aire maximale de la maille
125 aire_max = 0.
126 DO ij = 1, ip1jmp1
127 aire_max = max(aire_max, aire(ij))
128 END DO
129 ! factt = pas de temps en fraction de jour
130 factt = dtvr*iperiod/daysec
131
132 CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
133 CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
134 CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
135 CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
136 CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
137
138 CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
139 CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
140 CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
141
142 ! Cas ou on force exactement par les variables analysees
143 ELSE
144 alpha_t = 0.
145 alpha_u = 0.
146 alpha_v = 0.
147 alpha_p = 0.
148 ! physic=.false.
149 END IF
150
151 itau_test = 1001
152 step_rea = 1
153 count_no_rea = 0
154 ncidpl = -99
155
156 ! itau_test montre si l'importation a deja ete faite au rang itau
157 ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
158 if (guide_u) then
159 if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
160 endif
161
162 if (guide_v) then
163 if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
164 endif
165
166 if (guide_T) then
167 if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
168 endif
169
170 if (guide_Q) then
171 if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
172 endif
173
174 IF (ncep) THEN
175 status = nf90_inq_dimid(ncidpl, 'LEVEL', rid)
176 ELSE
177 status = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)
178 END IF
179 status = nf90_inquire_dimension(ncidpl, rid, len=nlev)
180 PRINT *, 'nlev', nlev
181 rcod = nf90_close(ncidpl)
182 ! Lecture du premier etat des reanalyses.
183 CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
184 masserea2, psrea2, 1, nlev)
185 qrea2(:, :) = max(qrea2(:, :), 0.1)
186
187 ! Debut de l'integration temporelle:
188 END IF ! first
189
190 ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
191
192 ditau = real(itau)
193 dday_step = real(day_step)
194 WRITE (*, *) 'ditau, dday_step'
195 WRITE (*, *) ditau, dday_step
196 toto = 4*ditau/dday_step
197 reste = toto - aint(toto)
198
199 IF (reste==0.) THEN
200 IF (itau_test==itau) THEN
201 WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
202 STOP
203 ELSE
204 vcovrea1(:, :) = vcovrea2(:, :)
205 ucovrea1(:, :) = ucovrea2(:, :)
206 tetarea1(:, :) = tetarea2(:, :)
207 qrea1(:, :) = qrea2(:, :)
208
209 PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
210 count_no_rea, ' non lectures'
211 step_rea = step_rea + 1
212 itau_test = itau
213 CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
214 qrea2, masserea2, psrea2, 1, nlev)
215 qrea2(:, :) = max(qrea2(:, :), 0.1)
216 factt = dtvr*iperiod/daysec
217 ztau(:) = factt/max(alpha_t(:), 1.E-10)
218 CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
219 CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
220 CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
221 CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
222 CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
223 CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
224 CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
225 CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
226 CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
227 CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
228 CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
229
230 CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')
231
232 END IF
233 ELSE
234 count_no_rea = count_no_rea + 1
235 END IF
236
237 ! Guidage
238 ! x_gcm = a * x_gcm + (1-a) * x_reanalyses
239
240 IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
241
242 ditau = real(itau)
243 dday_step = real(day_step)
244
245 tau = 4*ditau/dday_step
246 tau = tau - aint(tau)
247
248 ! ucov
249 IF (guide_u) THEN
250 DO l = 1, llm
251 DO ij = 1, ip1jmp1
252 a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)
253 ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a
254 IF (first .AND. ini_anal) ucov(ij, l) = a
255 END DO
256 END DO
257 END IF
258
259 IF (guide_t) THEN
260 DO l = 1, llm
261 DO ij = 1, ip1jmp1
262 a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)
263 teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a
264 IF (first .AND. ini_anal) teta(ij, l) = a
265 END DO
266 END DO
267 END IF
268
269 ! P
270 IF (guide_p) THEN
271 DO ij = 1, ip1jmp1
272 a = (1.-tau)*psrea1(ij) + tau*psrea2(ij)
273 ps(ij) = (1.-alpha_p(ij))*ps(ij) + alpha_p(ij)*a
274 IF (first .AND. ini_anal) ps(ij) = a
275 END DO
276 forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
277 CALL massdair(p, masse)
278 END IF
279
280 ! q
281 IF (guide_q) THEN
282 DO l = 1, llm
283 DO ij = 1, ip1jmp1
284 a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)
285 ! hum relative en % -> hum specif
286 a = qsat(ij, l)*a*0.01
287 q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a
288 IF (first .AND. ini_anal) q(ij, l) = a
289 END DO
290 END DO
291 END IF
292
293 ! vcov
294 IF (guide_v) THEN
295 DO l = 1, llm
296 DO ij = 1, ip1jm
297 a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)
298 vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a
299 IF (first .AND. ini_anal) vcov(ij, l) = a
300 END DO
301 IF (first .AND. ini_anal) vcov(ij, l) = a
302 END DO
303 END IF
304
305 first = .FALSE.
306
307 END SUBROUTINE guide
308
309 END MODULE guide_m

  ViewVC Help
Powered by ViewVC 1.1.21