Context Navigation

read_reanalyse.F @ 1261

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

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

Download in other formats:

Original Format