Context Navigation

read_reanalyse.F @ 1097

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

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

Download in other formats:

Original Format