Context Navigation

read_reanalyse.F @ 1065

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

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

Download in other formats:

Original Format