Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

read_reanalyse.F @ 1864

Last change on this file since 1864 was 703, checked in by (none), 18 years ago
This commit was manufactured by cvs2svn to create branch 'V3_test'.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 14.0 KB

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

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

Download in other formats:

Original Format