Context Navigation

read_reanalyse.F @ 1097

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

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

Download in other formats:

Original Format