Context Navigation

read_reanalyse.F @ 1086

Last change on this file since 1086 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.7 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
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,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*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 plnc(iip1,jjp1,nlevnc)
52	real*4 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
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	! endif (first)
131	endif
132	print*,'ok1'
133
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	#ifdef NC_DOUBLE
166	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unc)
167	#else
168	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
169	#endif
170	print*,'WARNING!!! Correction bidon pour palier a un '
171	print*,'probleme dans la creation des fichiers nc'
172	call correctbid(iim,jjp1*nlevnc,unc)
173	call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')
174	endif
175
176	c Temperature
177	c -----------
178
179	print*,'ncidt=',ncidt,'varidt=',varidt,'start=',start
180	print*,'count=',count
181	if (guide_T) then
182	#ifdef NC_DOUBLE
183	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnc)
184	#else
185	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc)
186	#endif
187	call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 AAA ')
188	call correctbid(iim,jjp1*nlevnc,tnc)
189	call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 BBB ')
190	endif
191
192	c Humidite
193	c --------
194
195	if (guide_Q) then
196	#ifdef NC_DOUBLE
197	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,Qnc)
198	#else
199	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)
200	#endif
201	call correctbid(iim,jjp1*nlevnc,Qnc)
202	call dump2d(iip1,jjp1,Qnc,'QNC COUCHE 1 ')
203	endif
204
205	count(2)=jjm
206	c Vent meridien
207	c -------------
208
209	if (guide_v) then
210	#ifdef NC_DOUBLE
211	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnc)
212	#else
213	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)
214	#endif
215	call correctbid(iim,jjm*nlevnc,vnc)
216	call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ')
217	endif
218
219	start(3)=timestep
220	start(4)=0
221	count(2)=jjp1
222	count(3)=1
223	count(4)=0
224
225	c Pression de surface
226	c -------------------
227	psnc(:,:)=0.
228	if ((guide_P).OR.(guide_modele)) then
229	#ifdef NC_DOUBLE
230	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc)
231	#else
232	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
233	#endif
234	call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')
235	call correctbid(iim,jjp1,psnc)
236	endif
237
238	c Calcul de la pression aux differents niveaux
239	c --------------------------------------------
240	if (guide_modele) then
241	#ifdef NC_DOUBLE
242	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
243	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
244	#else
245	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
246	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
247	#endif
248	else
249	#ifdef NC_DOUBLE
250	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
251	#else
252	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
253	#endif
254	c conversion en Pascals
255	apnc=apnc*100.
256	bpnc(:)=0.
257	endif
258	do i=1,iip1
259	do j=1,jjp1
260	do l=1,nlevnc
261	plnc(i,j,l)=apnc(l)+bpnc(l)*psnc(i,j)
262	enddo
263	enddo
264	enddo
265	if (first) then
266	print*,'Verification ordre niveaux verticaux'
267	print*,'LMDZ :'
268	do l=1,llm
269	print*,'PL(',l,')=',(ap(l)+ap(l+1))/2.
270	$ +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
271	enddo
272	print*,'Fichiers guidage'
273	do l=1,nlevnc
274	print*,'PL(',l,')=',plnc(1,1,l)
275	enddo
276	if (guide_u) then
277	do l=1,nlevnc
278	print*,'U(',l,')=',unc(1,1,l)
279	enddo
280	endif
281	if (guide_T) then
282	do l=1,nlevnc
283	print*,'T(',l,')=',tnc(1,1,l)
284	enddo
285	endif
286	print *,'inversion de l''ordre: ok_invertp=',ok_invertp
287	endif
288
289	c -----------------------------------------------------------------
290	c Interpollation verticale sur les niveaux modele
291	c -----------------------------------------------------------------
292	call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,plnc,u,v,t,Q
293	s ,ps,masse,pk)
294
295	call dump2d(iip1,jjm,v,'V COUCHE APRES ')
296
297
298	c -----------------------------------------------------------------
299	c Passage aux variables du modele (vents covariants, temperature
300	c potentielle et humidite specifique)
301	c -----------------------------------------------------------------
302	call nat2gcm(u,v,t,Q,pk,u,v,t,Q)
303	print*,'TIMESTEP ',timestep
304	if(mode.ne.1) stop'mode pas egal 0'
305	c call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ')
306
307	c Lignes introduites a une epoque pour un probleme oublie...
308	c do l=1,llm
309	c do i=1,iip1
310	c v(i,1,l)=0.
311	c v(i,jjm,l)=0.
312	c enddo
313	c enddo
314	first=.false.
315
316	return
317	end
318
319
320	c===========================================================================
321	subroutine reanalyse2nat(nlevnc,psi
322	s ,unc,vnc,tnc,qnc,psnc,plnc,u,v,t,q
323	s ,ps,masse,pk)
324	c===========================================================================
325
326	c -----------------------------------------------------------------
327	c Inversion Nord/sud de la grille + interpollation sur les niveaux
328	c verticaux du modele.
329	c -----------------------------------------------------------------
330
331	implicit none
332
333	#include "dimensions.h"
334	#include "paramet.h"
335	#include "comgeom2.h"
336	#include "comvert.h"
337	#include "comconst.h"
338	#include "guide.h"
339
340	integer nlevnc
341	real psi(iip1,jjp1)
342	real u(iip1,jjp1,llm),v(iip1,jjm,llm)
343	real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
344
345	real plnc(iip1,jjp1,nlevnc)
346	real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
347	real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
348	real qnc(iip1,jjp1,nlevnc)
349
350	real zu(iip1,jjp1,llm),zv(iip1,jjm,llm)
351	real zt(iip1,jjp1,llm),zq(iip1,jjp1,llm)
352
353	real pext(iip1,jjp1,llm)
354	real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm)
355	real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm)
356	real plsnc(iip1,jjp1,llm)
357
358	REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
359	real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1)
360	real pkf(iip1,jjp1,llm)
361	real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm)
362	real prefkap,unskap
363
364
365	integer i,j,l
366
367
368	c -----------------------------------------------------------------
369	c calcul de la pression au milieu des couches.
370	c -----------------------------------------------------------------
371
372	CALL pression( ip1jmp1, ap, bp, psi, p )
373	call massdair(p,masse)
374	CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
375
376	c .... Calcul de pls , pression au milieu des couches ,en Pascals
377	unskap=1./kappa
378	prefkap = preff ** kappa
379	c PRINT *,' Pref kappa unskap ',preff,kappa,unskap
380	DO l = 1, llm
381	DO j=1,jjp1
382	DO i =1, iip1
383	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
384	ENDDO
385	ENDDO
386	ENDDO
387
388
389	c -----------------------------------------------------------------
390	c calcul des pressions pour les grilles u et v
391	c -----------------------------------------------------------------
392
393	do l=1,llm
394	do j=1,jjp1
395	do i=1,iip1
396	pext(i,j,l)=pls(i,j,l)*aire(i,j)
397	enddo
398	enddo
399	enddo
400	call massbar(pext, pbarx, pbary )
401	do l=1,llm
402	do j=1,jjp1
403	do i=1,iip1
404	plunc(i,jjp1+1-j,l)=pbarx(i,j,l)/aireu(i,j)
405	plsnc(i,jjp1+1-j,l)=pls(i,j,l)
406	enddo
407	enddo
408	enddo
409	do l=1,llm
410	do j=1,jjm
411	do i=1,iip1
412	plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev(i,j)
413	enddo
414	enddo
415	enddo
416
417	c -----------------------------------------------------------------
418
419	if (guide_P) then
420	do j=1,jjp1
421	do i=1,iim
422	ps(i,j)=psnc(i,jjp1+1-j)
423	enddo
424	ps(iip1,j)=ps(1,j)
425	enddo
426	endif
427
428
429	c -----------------------------------------------------------------
430	call pres2lev(unc,zu,nlevnc,llm,plnc,plunc,iip1,jjp1,ok_invertp)
431	call pres2lev(vnc,zv,nlevnc,llm,plnc,plvnc,iip1,jjm,ok_invertp )
432	call pres2lev(tnc,zt,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)
433	call pres2lev(qnc,zq,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)
434
435	c call dump2d(iip1,jjp1,ps,'PS ')
436	c call dump2d(iip1,jjp1,psu,'PS ')
437	c call dump2d(iip1,jjm,psv,'PS ')
438	c Inversion Nord/Sud
439	do l=1,llm
440	do j=1,jjp1
441	do i=1,iim
442	u(i,j,l)=zu(i,jjp1+1-j,l)
443	t(i,j,l)=zt(i,jjp1+1-j,l)
444	q(i,j,l)=zq(i,jjp1+1-j,l)
445	enddo
446	u(iip1,j,l)=u(1,j,l)
447	t(iip1,j,l)=t(1,j,l)
448	q(iip1,j,l)=q(1,j,l)
449	enddo
450	enddo
451
452	do l=1,llm
453	do j=1,jjm
454	do i=1,iim
455	v(i,j,l)=zv(i,jjm+1-j,l)
456	enddo
457	v(iip1,j,l)=v(1,j,l)
458	enddo
459	enddo
460
461	return
462	end
463
464	c===========================================================================
465	subroutine nat2gcm(u,v,t,rh,pk,ucov,vcov,teta,q)
466	c===========================================================================
467
468	implicit none
469
470	#include "dimensions.h"
471	#include "paramet.h"
472	#include "comgeom2.h"
473	#include "comconst.h"
474	#include "comvert.h"
475	#include "guide.h"
476
477	real u(iip1,jjp1,llm),v(iip1,jjm,llm)
478	real t(iip1,jjp1,llm),pk(iip1,jjp1,llm),rh(iip1,jjp1,llm)
479	real ps(iip1,jjp1)
480
481	real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm)
482	real teta(iip1,jjp1,llm),q(iip1,jjp1,llm)
483
484	real pres(iip1,jjp1,llm),qsat(iip1,jjp1,llm)
485
486	real unskap
487
488	integer i,j,l
489
490
491	print*,'Entree dans nat2gcm'
492	c ucov(:,:,:)=0.
493	c do l=1,llm
494	c ucov(:,2:jjm,l)=u(:,2:jjm,l)*cu(:,2:jjm)
495	c enddo
496	c ucov(iip1,:,:)=ucov(1,:,:)
497
498	c teta(:,:,:)=t(:,:,:)*cpp/pk(:,:,:)
499	c teta(iip1,:,:)=teta(1,:,:)
500
501	c calcul de ucov et de la temperature potentielle
502	do l=1,llm
503	do j=1,jjp1
504	do i=1,iim
505	ucov(i,j,l)=u(i,j,l)*cu(i,j)
506	teta(i,j,l)=t(i,j,l)*cpp/pk(i,j,l)
507	enddo
508	ucov(iip1,j,l)=ucov(1,j,l)
509	teta(iip1,j,l)=teta(1,j,l)
510	enddo
511	do i=1,iip1
512	ucov(i,1,l)=0.
513	ucov(i,jjp1,l)=0.
514	teta(i,1,l)=teta(1,1,l)
515	teta(i,jjp1,l)=teta(1,jjp1,l)
516	enddo
517	enddo
518
519	c calcul de ucov
520	do l=1,llm
521	do j=1,jjm
522	do i=1,iim
523	vcov(i,j,l)=v(i,j,l)*cv(i,j)
524	enddo
525	vcov(iip1,j,l)=vcov(1,j,l)
526	enddo
527	enddo
528
529	c call dump2d(iip1,jjp1,teta,'TETA EN BAS ')
530	c call dump2d(iip1,jjp1,teta(1,1,llm),'TETA EN HAUT ')
531
532	c Humidite relative -> specifique
533	c -------------------------------
534	if (guide_hr) then
535	c FINALEMENT ON GUIDE EN HUMIDITE RELATIVE
536	print*,'calcul de unskap'
537	unskap = 1./ kappa
538	print*,'calcul de pres'
539	pres(:,:,:)=preff(pk(:,:,:)/cpp)*unskap
540	print*,'calcul de qsat'
541	call q_sat(iip1jjp1llm,t,pres,qsat)
542	print*,'calcul de q'
543	c ATTENTION : humidites relatives en %
544	rh(:,:,:)=max(rh(:,:,:)*0.01,1.e-6)
545	q(:,:,:)=qsat(:,:,:)*rh(:,:,:)
546	print*,'calcul de q OK'
547	call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1 ')
548	call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE ')
549	else
550	q(:,:,:)=rh(:,:,:)
551	print*,'calcul de q OK'
552	call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1 ')
553	endif
554
555	return
556	end
557
558
559
560	c===========================================================================
561	subroutine correctbid(iim,nl,x)
562	c===========================================================================
563	integer iim,nl
564	real x(iim+1,nl)
565	integer i,l
566	real zz
567
568	do l=1,nl
569	do i=2,iim-1
570	if(abs(x(i,l)).gt.1.e10) then
571	zz=0.5*(x(i-1,l)+x(i+1,l))
572	c print*,'correction ',i,l,x(i,l),zz
573	x(i,l)=zz
574	endif
575	enddo
576	enddo
577	return
578	end

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

Download in other formats:

Original Format