/[lmdze]/trunk/libf/dyn3d/read_reanalyse.f
ViewVC logotype

Contents of /trunk/libf/dyn3d/read_reanalyse.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 20 - (show annotations)
Wed Oct 15 16:19:57 2008 UTC (15 years, 7 months ago) by guez
File size: 13994 byte(s)
Deleted argument "presnivs" of "physiq", "ini_histhf", "ini_histhf3d",
"ini_histday", "ini_histins", "ini_histrac", "phytrac". Access it from
"comvert" instead.

Replaced calls to NetCDF Fortran 77 interface by calls to Fortran 90
interface or to NetCDF95.

Procedure "gr_phy_write_3d" now works with a variable of arbitrary
size in the second dimension.

Annotated use statements with "only" clause.

Replaced calls to NetCDF interface version 2 by calls to Fortran 90
interface in "guide.f90" and "read_reanalyse.f".

In "write_histrac", replaced calls to "gr_fi_ecrit" by calls to
"gr_phy_write_2d" and "gr_phy_write_3d".

1 !
2 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/read_reanalyse.F,v 1.3 2005/04/15 12:31:21 lmdzadmin Exp $
3 !
4 c
5 c
6 subroutine read_reanalyse(timestep,psi
7 s ,u,v,t,q,masse,ps,mode,nlevnc)
8
9 c mode=0 variables naturelles
10 c mode=1 variabels GCM
11
12 c -----------------------------------------------------------------
13 c Declarations
14 c -----------------------------------------------------------------
15 use dimens_m
16 use paramet_m
17 use comvert
18 use comgeom
19 use guide_m
20 use netcdf
21
22 IMPLICIT NONE
23
24 c common
25 c ------
26
27 include "netcdf.inc"
28
29
30 c arguments
31 c ---------
32 integer nlevnc
33 integer timestep,mode,l
34
35 real psi(iip1,jjp1)
36 real u(iip1,jjp1,llm),v(iip1,jjm,llm)
37 real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
38 real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
39
40
41 c local
42 c -----
43 integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
44 integer ncidpl
45 integer varidpl,ncidQ,varidQ
46 save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
47 save ncidpl
48 save varidpl,ncidQ,varidQ
49
50 real*4 unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
51 real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
52 real*4 Qnc(iip1,jjp1,nlevnc)
53 real*4 pl(nlevnc)
54
55 integer start(4),count(4),status
56
57 real rcode
58 logical first
59 save first
60
61 data first/.true./
62
63
64
65 c -----------------------------------------------------------------
66 c Initialisation de la lecture des fichiers
67 c -----------------------------------------------------------------
68 if (first) then
69 ncidpl=-99
70 print*,'Intitialisation de read reanalsye'
71
72 c Vent zonal
73 if (guide_u) then
74 rcode=nf90_open('u.nc',nf90_nowrite,ncidu)
75 rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
76 print*,'ncidu,varidu',ncidu,varidu
77 if (ncidpl.eq.-99) ncidpl=ncidu
78 endif
79
80 c Vent meridien
81 if (guide_v) then
82 rcode=nf90_open('v.nc',nf90_nowrite,ncidv)
83 rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
84 print*,'ncidv,varidv',ncidv,varidv
85 if (ncidpl.eq.-99) ncidpl=ncidv
86 endif
87
88 c Temperature
89 if (guide_T) then
90 rcode=nf90_open('T.nc',nf90_nowrite,ncidt)
91 rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
92 print*,'ncidt,varidt',ncidt,varidt
93 if (ncidpl.eq.-99) ncidpl=ncidt
94 endif
95
96 c Humidite
97 if (guide_Q) then
98 rcode=nf90_open('hur.nc',nf90_nowrite,ncidQ)
99 rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
100 print*,'ncidQ,varidQ',ncidQ,varidQ
101 if (ncidpl.eq.-99) ncidpl=ncidQ
102 endif
103
104 c Pression de surface
105 if (guide_P) then
106 rcode=nf90_open('ps.nc',nf90_nowrite,ncidps)
107 rcode = nf90_inq_varid(ncidps, 'SP', varidps)
108 print*,'ncidps,varidps',ncidps,varidps
109 endif
110
111 c Coordonnee verticale
112 if (ncep) then
113 print*,'Vous etes entrain de lire des donnees NCEP'
114 rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
115 else
116 print*,'Vous etes entrain de lire des donnees ECMWF'
117 rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
118 endif
119 print*,'ncidu,varidpl',ncidu,varidpl
120 endif
121 print*,'ok1'
122
123 c Niveaux de pression
124 print*,'WARNING!!! Il n y a pas de test de coherence'
125 print*,'sur le nombre de niveaux verticaux dans le fichier nc'
126 status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,pl)
127 c passage en pascal
128 pl(:)=100.*pl(:)
129 if (first) then
130 do l=1,nlevnc
131 print*,'PL(',l,')=',pl(l)
132 enddo
133 endif
134
135 c -----------------------------------------------------------------
136 c lecture des champs u, v, T, ps
137 c -----------------------------------------------------------------
138
139 c dimensions pour les champs scalaires et le vent zonal
140 c -----------------------------------------------------
141
142 start(1)=1
143 start(2)=1
144 start(3)=1
145 start(4)=timestep
146
147 count(1)=iip1
148 count(2)=jjp1
149 count(3)=nlevnc
150 count(4)=1
151
152 c mise a zero des tableaux
153 c ------------------------
154 unc(:,:,:)=0.
155 vnc(:,:,:)=0.
156 tnc(:,:,:)=0.
157 Qnc(:,:,:)=0.
158
159 c Vent zonal
160 c ----------
161
162 if (guide_u) then
163 print*,'avant la lecture de UNCEP nd de niv:',nlevnc
164 status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
165 c call dump2d(iip1,jjp1,unc,'VENT NCEP ')
166 c call dump2d(iip1,40,unc(1,1,nlevnc),'VENT NCEP ')
167 print*,'WARNING!!! Correction bidon pour palier a un '
168 print*,'probleme dans la creation des fichiers nc'
169 call correctbid(iim,jjp1*nlevnc,unc)
170 call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')
171 endif
172
173 c Temperature
174 c -----------
175
176 print*,'ncidt=',ncidt,'varidt=',varidt,'start=',start
177 print*,'count=',count
178 if (guide_T) then
179 status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc)
180 call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 AAA ')
181 call correctbid(iim,jjp1*nlevnc,tnc)
182 call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 BBB ')
183 endif
184
185 c Humidite
186 c --------
187
188 if (guide_Q) then
189 status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)
190 call correctbid(iim,jjp1*nlevnc,Qnc)
191 call dump2d(iip1,jjp1,Qnc,'QNC COUCHE 1 ')
192 endif
193
194 count(2)=jjm
195 c Vent meridien
196 c -------------
197
198 if (guide_v) then
199 status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)
200 call correctbid(iim,jjm*nlevnc,vnc)
201 call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ')
202 endif
203
204 start(3)=timestep
205 start(4)=0
206 count(2)=jjp1
207 count(3)=1
208 count(4)=0
209
210 c Pression de surface
211 c -------------------
212
213 if (guide_P) then
214 status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
215 call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')
216 call correctbid(iim,jjp1,psnc)
217 endif
218
219
220
221 c -----------------------------------------------------------------
222 c Interpollation verticale sur les niveaux modele
223 c -----------------------------------------------------------------
224 call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,pl,u,v,t,Q
225 s ,ps,masse,pk)
226
227 call dump2d(iip1,jjm,v,'V COUCHE APRES ')
228
229
230 c -----------------------------------------------------------------
231 c Passage aux variables du modele (vents covariants, temperature
232 c potentielle et humidite specifique)
233 c -----------------------------------------------------------------
234 call nat2gcm(u,v,t,Q,pk,u,v,t,Q)
235 print*,'TIMESTEP ',timestep
236 if(mode.ne.1) stop'mode pas egal 0'
237 c call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ')
238
239 c Lignes introduites a une epoque pour un probleme oublie...
240 c do l=1,llm
241 c do i=1,iip1
242 c v(i,1,l)=0.
243 c v(i,jjm,l)=0.
244 c enddo
245 c enddo
246 first=.false.
247
248 return
249 end
250
251
252 c===========================================================================
253 subroutine reanalyse2nat(nlevnc,psi
254 s ,unc,vnc,tnc,qnc,psnc,pl,u,v,t,q
255 s ,ps,masse,pk)
256 c===========================================================================
257
258 c -----------------------------------------------------------------
259 c Inversion Nord/sud de la grille + interpollation sur les niveaux
260 c verticaux du modele.
261 c -----------------------------------------------------------------
262
263 use dimens_m
264 use paramet_m
265 use comconst
266 use comvert
267 use comgeom
268 use exner_hyb_m, only: exner_hyb
269 use guide_m
270 use pression_m, only: pression
271
272 implicit none
273
274
275 integer nlevnc
276 real psi(iip1,jjp1)
277 real u(iip1,jjp1,llm),v(iip1,jjm,llm)
278 real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
279
280 real pl(nlevnc)
281 real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
282 real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
283 real qnc(iip1,jjp1,nlevnc)
284
285 real zu(iip1,jjp1,llm),zv(iip1,jjm,llm)
286 real zt(iip1,jjp1,llm),zq(iip1,jjp1,llm)
287
288 real pext(iip1,jjp1,llm)
289 real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm)
290 real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm)
291 real plsnc(iip1,jjp1,llm)
292
293 real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1)
294 real pkf(iip1,jjp1,llm)
295 real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm)
296 real prefkap,unskap
297
298
299 integer i,j,l
300
301
302 c -----------------------------------------------------------------
303 c calcul de la pression au milieu des couches.
304 c -----------------------------------------------------------------
305
306 CALL pression( ip1jmp1, ap, bp, psi, p )
307 call massdair(p,masse)
308 CALL exner_hyb(psi,p,pks,pk,pkf)
309
310 c .... Calcul de pls , pression au milieu des couches ,en Pascals
311 unskap=1./kappa
312 prefkap = preff ** kappa
313 c PRINT *,' Pref kappa unskap ',preff,kappa,unskap
314 DO l = 1, llm
315 DO j=1,jjp1
316 DO i =1, iip1
317 pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
318 ENDDO
319 ENDDO
320 ENDDO
321
322
323 c -----------------------------------------------------------------
324 c calcul des pressions pour les grilles u et v
325 c -----------------------------------------------------------------
326
327 do l=1,llm
328 do j=1,jjp1
329 do i=1,iip1
330 pext(i,j,l)=pls(i,j,l)*aire_2d(i,j)
331 enddo
332 enddo
333 enddo
334 call massbar(pext, pbarx, pbary )
335 do l=1,llm
336 do j=1,jjp1
337 do i=1,iip1
338 plunc(i,jjp1+1-j,l)=pbarx(i,j,l)/aireu_2d(i,j)
339 plsnc(i,jjp1+1-j,l)=pls(i,j,l)
340 enddo
341 enddo
342 enddo
343 do l=1,llm
344 do j=1,jjm
345 do i=1,iip1
346 plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev_2d(i,j)
347 enddo
348 enddo
349 enddo
350
351 c -----------------------------------------------------------------
352
353 if (guide_P) then
354 do j=1,jjp1
355 do i=1,iim
356 ps(i,j)=psnc(i,jjp1+1-j)
357 enddo
358 ps(iip1,j)=ps(1,j)
359 enddo
360 endif
361
362
363 c -----------------------------------------------------------------
364 call pres2lev(unc,zu,nlevnc,llm,pl,plunc,iip1,jjp1)
365 call pres2lev(vnc,zv,nlevnc,llm,pl,plvnc,iip1,jjm )
366 call pres2lev(tnc,zt,nlevnc,llm,pl,plsnc,iip1,jjp1)
367 call pres2lev(qnc,zq,nlevnc,llm,pl,plsnc,iip1,jjp1)
368
369 c call dump2d(iip1,jjp1,ps,'PS ')
370 c call dump2d(iip1,jjp1,psu,'PS ')
371 c call dump2d(iip1,jjm,psv,'PS ')
372 c Inversion Nord/Sud
373 do l=1,llm
374 do j=1,jjp1
375 do i=1,iim
376 u(i,j,l)=zu(i,jjp1+1-j,l)
377 t(i,j,l)=zt(i,jjp1+1-j,l)
378 q(i,j,l)=zq(i,jjp1+1-j,l)
379 enddo
380 u(iip1,j,l)=u(1,j,l)
381 t(iip1,j,l)=t(1,j,l)
382 q(iip1,j,l)=q(1,j,l)
383 enddo
384 enddo
385
386 do l=1,llm
387 do j=1,jjm
388 do i=1,iim
389 v(i,j,l)=zv(i,jjm+1-j,l)
390 enddo
391 v(iip1,j,l)=v(1,j,l)
392 enddo
393 enddo
394
395 return
396 end
397
398 c===========================================================================
399 subroutine nat2gcm(u,v,t,rh,pk,ucov,vcov,teta,q)
400 c===========================================================================
401
402 use dimens_m
403 use paramet_m
404 use comconst
405 use comvert
406 use comgeom
407 use q_sat_m, only: q_sat
408 use guide_m
409 implicit none
410
411
412 real u(iip1,jjp1,llm),v(iip1,jjm,llm)
413 real t(iip1,jjp1,llm),pk(iip1,jjp1,llm),rh(iip1,jjp1,llm)
414 real ps(iip1,jjp1)
415
416 real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm)
417 real teta(iip1,jjp1,llm),q(iip1,jjp1,llm)
418
419 real pres(iip1,jjp1,llm),qsat(iip1,jjp1,llm)
420
421 real unskap
422
423 integer i,j,l
424
425
426 print*,'Entree dans nat2gcm'
427 c ucov(:,:,:)=0.
428 c do l=1,llm
429 c ucov(:,2:jjm,l)=u(:,2:jjm,l)*cu_2d(:,2:jjm)
430 c enddo
431 c ucov(iip1,:,:)=ucov(1,:,:)
432
433 c teta(:,:,:)=t(:,:,:)*cpp/pk(:,:,:)
434 c teta(iip1,:,:)=teta(1,:,:)
435
436 c calcul de ucov et de la temperature potentielle
437 do l=1,llm
438 do j=1,jjp1
439 do i=1,iim
440 ucov(i,j,l)=u(i,j,l)*cu_2d(i,j)
441 teta(i,j,l)=t(i,j,l)*cpp/pk(i,j,l)
442 enddo
443 ucov(iip1,j,l)=ucov(1,j,l)
444 teta(iip1,j,l)=teta(1,j,l)
445 enddo
446 do i=1,iip1
447 ucov(i,1,l)=0.
448 ucov(i,jjp1,l)=0.
449 teta(i,1,l)=teta(1,1,l)
450 teta(i,jjp1,l)=teta(1,jjp1,l)
451 enddo
452 enddo
453
454 c calcul de ucov
455 do l=1,llm
456 do j=1,jjm
457 do i=1,iim
458 vcov(i,j,l)=v(i,j,l)*cv_2d(i,j)
459 enddo
460 vcov(iip1,j,l)=vcov(1,j,l)
461 enddo
462 enddo
463
464 c call dump2d(iip1,jjp1,teta,'TETA EN BAS ')
465 c call dump2d(iip1,jjp1,teta(1,1,llm),'TETA EN HAUT ')
466
467 c Humidite relative -> specifique
468 c -------------------------------
469 if (1.eq.0) then
470 c FINALEMENT ON GUIDE EN HUMIDITE RELATIVE
471 print*,'calcul de unskap'
472 unskap = 1./ kappa
473 print*,'calcul de pres'
474 pres(:,:,:)=preff*(pk(:,:,:)/cpp)**unskap
475 print*,'calcul de qsat'
476 qsat = q_sat(t, pres)
477 print*,'calcul de q'
478 c ATTENTION : humidites relatives en %
479 rh(:,:,:)=max(rh(:,:,:)*0.01,1.e-6)
480 q(:,:,:)=qsat(:,:,:)*rh(:,:,:)
481 print*,'calcul de q OK'
482
483 call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE ')
484 call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1 ')
485 endif
486
487
488 return
489 end
490
491
492
493 c===========================================================================
494 subroutine correctbid(iim,nl,x)
495 c===========================================================================
496 integer iim,nl
497 real x(iim+1,nl)
498 integer i,l
499 real zz
500
501 do l=1,nl
502 do i=2,iim-1
503 if(abs(x(i,l)).gt.1.e10) then
504 zz=0.5*(x(i-1,l)+x(i+1,l))
505 c print*,'correction ',i,l,x(i,l),zz
506 x(i,l)=zz
507 endif
508 enddo
509 enddo
510 return
511 end

  ViewVC Help
Powered by ViewVC 1.1.21