Context Navigation

read_reanalyse.F @ 1114

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

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

Download in other formats:

Original Format