Context Navigation

read_reanalyse.F @ 617

Last change on this file since 617 was 617, checked in by lmdzadmin, 19 years ago

corrections pour le mode guide :

mode non zoome
guidage de certains champs seulement

MAF

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 14.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format