Context Navigation

read_reanalyse.F @ 3947

Last change on this file since 3947 was 753, checked in by lsce, 18 years ago
ACo : correction pour le mode guide en parallele
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 14.3 KB

Rev	Line
[630]	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
[753]	12	USE parallel
[630]	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
[753]	45	save ncidpl,ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
[630]	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
[753]	58	INTEGER ierr
[630]	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
[753]	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
[630]	77	endif
	78
	79	c Vent meridien
	80	if (guide_v) then
[753]	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
[630]	85	endif
[753]	86
[630]	87
	88	c Temperature
	89	if (guide_T) then
[753]	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
[630]	94	endif
[753]	95
[630]	96	c Humidite
	97	if (guide_Q) then
[753]	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
[630]	102	endif
[753]	103
[630]	104	c Pression de surface
	105	if (guide_P) then
[753]	106	ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
	107	varidps=NCVID(ncidps,'SP',rcode)
	108	print*,'ncidps,varidps',ncidps,varidps
[630]	109	endif
[753]	110
[630]	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
[753]	139
[630]	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
[753]	166
	167	c mise a zero des tableaux
	168	c ------------------------
	169	unc(:,:,:)=0.
	170	vnc(:,:,:)=0.
	171	tnc(:,:,:)=0.
	172	Qnc(:,:,:)=0.
	173
[630]	174	c Vent zonal
	175	c ----------
	176
	177	if (guide_u) then
[753]	178	print*,'avant la lecture de UNCEP nd de niv:',nlevnc
	179
[630]	180	#ifdef NC_DOUBLE
[753]	181	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unc)
[630]	182	#else
[753]	183	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
[630]	184	#endif
	185	c call dump2d(iip1,jjp1,unc,'VENT NCEP ')
	186	c call dump2d(iip1,40,unc(1,1,nlevnc),'VENT NCEP ')
[753]	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 ')
[630]	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
[753]	200	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnc)
[630]	201	#else
[753]	202	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc)
[630]	203	#endif
[753]	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 ')
[630]	207	endif
	208
	209	c Humidite
	210	c --------
	211
	212	if (guide_Q) then
	213	#ifdef NC_DOUBLE
[753]	214	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,Qnc)
[630]	215	#else
[753]	216	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)
[630]	217	#endif
[753]	218	call correctbid(iim,jjp1*nlevnc,Qnc)
	219	call dump2d(iip1,jjp1,Qnc,'QNC COUCHE 1 ')
[630]	220	endif
	221
	222	count(2)=jjm
	223	c Vent meridien
	224	c -------------
	225
	226	if (guide_v) then
	227	#ifdef NC_DOUBLE
[753]	228	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnc)
[630]	229	#else
[753]	230	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)
[630]	231	#endif
[753]	232	call correctbid(iim,jjm*nlevnc,vnc)
	233	call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ')
[630]	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
[753]	247	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc)
[630]	248	#else
[753]	249	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
[630]	250	#endif
[753]	251	call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')
	252	call correctbid(iim,jjp1,psnc)
[630]	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