/[lmdze]/trunk/Sources/dyn3d/Guide/Read_reanalyse/read_reanalyse.f
ViewVC logotype

Annotation of /trunk/Sources/dyn3d/Guide/Read_reanalyse/read_reanalyse.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21