Context Navigation

read_reanalyse.F @ 701

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

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

Download in other formats:

Original Format