Context Navigation

read_reanalyse.F @ 1198

Last change on this file since 1198 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: 15.9 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 -----------------------------------------------------------------
[1146]	15	use netcdf
	16
[524]	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
	30	c arguments
	31	c ---------
	32	integer nlevnc
[1046]	33	integer timestep,mode
[524]	34
	35	real psi(iip1,jjp1)
	36	real u(iip1,jjp1,llm),v(iip1,jjm,llm)
	37	real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
	38	real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
	39
	40
	41	c local
	42	c -----
	43	integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
	44	integer ncidpl
[1046]	45	integer varidpl,ncidQ,varidQ,varidap,varidbp
[524]	46	save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
[1046]	47	save ncidpl,varidap,varidbp
[524]	48	save varidpl,ncidQ,varidQ
	49
[1122]	50	real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
	51	real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
	52	real Qnc(iip1,jjp1,nlevnc)
	53	real plnc(iip1,jjp1,nlevnc)
	54	real apnc(nlevnc),bpnc(nlevnc)
[524]	55
	56	integer start(4),count(4),status
[1046]	57	integer i,j,l
[524]	58
	59	real rcode
	60	logical first
	61	save first
	62
	63	data first/.true./
	64
	65	c -----------------------------------------------------------------
	66	c Initialisation de la lecture des fichiers
	67	c -----------------------------------------------------------------
	68	if (first) then
	69	ncidpl=-99
	70	print*,'Intitialisation de read reanalsye'
	71
[1046]	72	c Niveaux de pression si non constants
	73	if (guide_modele) then
	74	print *,'Vous êtes entrain de lire des données sur
	75	. niveaux modèle'
[1146]	76	rcode=nf90_open('apbp.nc',nf90_nowrite,ncidpl)
	77	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
	78	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
[1046]	79	print*,'ncidpl,varidap',ncidpl,varidap
	80	endif
	81
[524]	82	c Vent zonal
	83	if (guide_u) then
[1146]	84	rcode=nf90_open('u.nc',nf90_nowrite,ncidu)
	85	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
	86	print*,'ncidu,varidu',ncidu,varidu
	87	if (ncidpl.eq.-99) ncidpl=ncidu
[524]	88	endif
	89
	90	c Vent meridien
	91	if (guide_v) then
[1146]	92	rcode=nf90_open('v.nc',nf90_nowrite,ncidv)
	93	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
[524]	94	print*,'ncidv,varidv',ncidv,varidv
[617]	95	if (ncidpl.eq.-99) ncidpl=ncidv
[524]	96	endif
	97
	98	c Temperature
	99	if (guide_T) then
[1146]	100	rcode=nf90_open('T.nc',nf90_nowrite,ncidt)
	101	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
[524]	102	print*,'ncidt,varidt',ncidt,varidt
[617]	103	if (ncidpl.eq.-99) ncidpl=ncidt
[524]	104	endif
	105
	106	c Humidite
	107	if (guide_Q) then
[1146]	108	rcode=nf90_open('hur.nc',nf90_nowrite,ncidQ)
	109	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
[524]	110	print*,'ncidQ,varidQ',ncidQ,varidQ
[617]	111	if (ncidpl.eq.-99) ncidpl=ncidQ
[524]	112	endif
	113
	114	c Pression de surface
[1046]	115	if ((guide_P).OR.(guide_modele)) then
[1146]	116	rcode=nf90_open('ps.nc',nf90_nowrite,ncidps)
	117	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
[524]	118	print*,'ncidps,varidps',ncidps,varidps
	119	endif
	120
[618]	121	c Coordonnee verticale
[1046]	122	if (.not.guide_modele) then
[1146]	123	if (ncep) then
	124	print*,'Vous etes entrain de lire des donnees NCEP'
	125	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
	126	else
	127	print*,'Vous etes entrain de lire des donnees ECMWF'
	128	rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
	129	endif
	130	print*,'ncidpl,varidpl',ncidpl,varidpl
[524]	131	endif
[1046]	132	! endif (first)
[524]	133	endif
	134	print*,'ok1'
	135
	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.
[1046]	160	plnc(:,:,:)=0.
[617]	161
[524]	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	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 ')
	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
	185	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnc)
	186	#else
	187	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc)
	188	#endif
	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 ')
	192	endif
	193
	194	c Humidite
	195	c --------
	196
	197	if (guide_Q) then
	198	#ifdef NC_DOUBLE
	199	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,Qnc)
	200	#else
	201	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)
	202	#endif
	203	call correctbid(iim,jjp1*nlevnc,Qnc)
	204	call dump2d(iip1,jjp1,Qnc,'QNC COUCHE 1 ')
	205	endif
	206
	207	count(2)=jjm
	208	c Vent meridien
	209	c -------------
	210
	211	if (guide_v) then
	212	#ifdef NC_DOUBLE
	213	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnc)
	214	#else
	215	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)
	216	#endif
	217	call correctbid(iim,jjm*nlevnc,vnc)
	218	call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ')
	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
[524]	231	#ifdef NC_DOUBLE
	232	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc)
	233	#else
	234	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
	235	#endif
	236	call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')
	237	call correctbid(iim,jjp1,psnc)
	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
[524]	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
[524]	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
[524]	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)
[524]	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)
[524]	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
[524]	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 ')
[524]	550	call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE ')
[1046]	551	else
	552	q(:,:,:)=rh(:,:,:)
	553	print*,'calcul de q OK'
[524]	554	call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1 ')
[1046]	555	endif
[524]	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