Context Navigation

read_reanalyse.F @ 1275

Last change on this file since 1275 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.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	USE parallel
13	use netcdf
14	c -----------------------------------------------------------------
15	c Declarations
16	c -----------------------------------------------------------------
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	c arguments
30	c ---------
31	integer nlevnc
32	integer timestep,mode
33
34	real psi(iip1,jjp1)
35	real u(iip1,jjp1,llm),v(iip1,jjm,llm)
36	real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
37	real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
38
39	c local
40	c -----
41	integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
42	integer ncidpl
43	integer varidpl,ncidQ,varidQ,varidap,varidbp
44	save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
45	save ncidpl,varidap,varidbp
46	save varidpl,ncidQ,varidQ
47
48	real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
49	real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
50	real Qnc(iip1,jjp1,nlevnc)
51	real plnc(iip1,jjp1,nlevnc)
52	real apnc(nlevnc),bpnc(nlevnc)
53
54	integer start(4),count(4),status
55	integer i,j,l
56
57	real rcode
58	logical first
59	save first
60	INTEGER ierr
61
62	data first/.true./
63
64	c -----------------------------------------------------------------
65	c Initialisation de la lecture des fichiers
66	c -----------------------------------------------------------------
67	if (first) then
68	ncidpl=-99
69	print*,'Intitialisation de read reanalsye'
70
71	c Niveaux de pression si non constants
72	if (guide_modele) then
73	print *,'Vous êtes entrain de lire des données sur
74	. niveaux modèle'
75	rcode=nf90_open('apbp.nc',nf90_nowrite,ncidpl)
76	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
77	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
78	print*,'ncidpl,varidap',ncidpl,varidap
79	endif
80
81	c Vent zonal
82	if (guide_u) then
83	rcode=nf90_open('u.nc',nf90_nowrite,ncidu)
84	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
85	print*,'ncidu,varidu',ncidu,varidu
86	if (ncidpl.eq.-99) ncidpl=ncidu
87	endif
88
89	c Vent meridien
90	if (guide_v) then
91	rcode=nf90_open('v.nc',nf90_nowrite,ncidv)
92	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
93	print*,'ncidv,varidv',ncidv,varidv
94	if (ncidpl.eq.-99) ncidpl=ncidv
95	endif
96
97	c Temperature
98	if (guide_T) then
99	rcode=nf90_open('T.nc',nf90_nowrite,ncidt)
100	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
101	print*,'ncidt,varidt',ncidt,varidt
102	if (ncidpl.eq.-99) ncidpl=ncidt
103	endif
104
105	c Humidite
106	if (guide_Q) then
107	rcode=nf90_open('hur.nc',nf90_nowrite,ncidQ)
108	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
109	print*,'ncidQ,varidQ',ncidQ,varidQ
110	if (ncidpl.eq.-99) ncidpl=ncidQ
111	endif
112
113	c Pression de surface
114	if ((guide_P).OR.(guide_modele)) then
115	rcode=nf90_open('ps.nc',nf90_nowrite,ncidps)
116	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
117	print*,'ncidps,varidps',ncidps,varidps
118	endif
119
120	c Coordonnee verticale
121	if (.not.guide_modele) then
122	if (ncep) then
123	print*,'Vous etes entrain de lire des donnees NCEP'
124	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
125	else
126	print*,'Vous etes entrain de lire des donnees ECMWF'
127	rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
128	endif
129	print*,'ncidpl,varidpl',ncidpl,varidpl
130	endif
131	print*,'ncidu,varidpl',ncidu,varidpl
132
133	endif
134	print*,'ok1'
135
136	c -----------------------------------------------------------------
137	c lecture des champs u, v, T, ps
138	c -----------------------------------------------------------------
139
140	c dimensions pour les champs scalaires et le vent zonal
141	c -----------------------------------------------------
142
143	start(1)=1
144	start(2)=1
145	start(3)=1
146	start(4)=timestep
147
148	count(1)=iip1
149	count(2)=jjp1
150	count(3)=nlevnc
151	count(4)=1
152
153	c mise a zero des tableaux
154	c ------------------------
155	unc(:,:,:)=0.
156	vnc(:,:,:)=0.
157	tnc(:,:,:)=0.
158	Qnc(:,:,:)=0.
159	plnc(:,:,:)=0.
160
161	c Vent zonal
162	c ----------
163
164	if (guide_u) then
165	print*,'avant la lecture de UNCEP nd de niv:',nlevnc
166
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