Context Navigation

read_reanalyse.F @ 1198

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

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

Download in other formats:

Original Format