/[lmdze]/trunk/dyn3d/guide.f
ViewVC logotype

Contents of /trunk/dyn3d/guide.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 88 - (show annotations)
Tue Mar 11 15:09:02 2014 UTC (10 years, 2 months ago) by guez
File size: 9821 byte(s)
Removed useless argument mode of subroutine read_reanalyse.

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

  ViewVC Help
Powered by ViewVC 1.1.21