Context Navigation

read_reanalyse.F @ 1171

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

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

Download in other formats:

Original Format