Context Navigation

limit_netcdf.F @ 5306

Last change on this file since 5306 was 425, checked in by lmdz, 22 years ago

Dans le cas "oldice" la routine "sea_ice" utilise seulement l'interpol.
grille_m Levan

LF

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 28.9 KB

Rev	Line
[259]	1	C
	2	C $Header$
	3	C
	4	SUBROUTINE limit_netcdf ( interbar, extrap, oldice )
	5	c
	6	IMPLICIT none
	7	c
	8	c-------------------------------------------------------------
	9	C Author : L. Fairhead
	10	C Date : 27/01/94
	11	C Objet : Construction des fichiers de conditions aux limites
	12	C pour le nouveau
	13	C modele a partir de fichiers de climatologie. Les deux
	14	C grilles doivent etre regulieres
	15	c
	16	c Modifie par z.x.li (le23mars1994)
	17	c Modifie par L. Fairhead (fairhead@lmd.jussieu.fr) septembre 1999
	18	c pour lecture netcdf dans LMDZ.3.3
	19	c Modifie par P;Le Van , juillet 2001
	20	c-------------------------------------------------------------
	21	c
	22	#include "dimensions.h"
	23	#include "paramet.h"
	24	#include "control.h"
	25	#include "logic.h"
	26	#include "netcdf.inc"
	27	#include "comvert.h"
	28	#include "comgeom2.h"
	29	#include "comconst.h"
	30	c
	31	c-----------------------------------------------------------------------
	32	LOGICAL interbar, extrap, oldice
	33
	34	INTEGER KIDIA, KFDIA, KLON, KLEV
	35	PARAMETER (KIDIA=1,KFDIA=iim*(jjm-1)+2,
	36	. KLON=KFDIA-KIDIA+1,KLEV=llm)
	37	c-----------------------------------------------------------------------
	38	REAL phy_nat(klon,360), phy_nat0(klon)
	39	REAL phy_alb(klon,360)
	40	REAL phy_sst(klon,360)
	41	REAL phy_bil(klon,360)
	42	REAL phy_rug(klon,360)
	43	REAL phy_ice(klon,360)
	44	c
	45	REAL masque(iip1,jjp1)
	46	REAL mask(iim,jjp1)
	47
	48	C Declarations pour le champ de depart
	49	INTEGER imdep, jmdep,lmdep
	50	INTEGER tbid
	51	PARAMETER ( tbid = 60 ) ! >52 semaines
	52	REAL timecoord(tbid)
	53	c
	54	REAL , ALLOCATABLE :: dlon_msk(:), dlat_msk(:)
	55	REAL , ALLOCATABLE :: lonmsk_ini(:), latmsk_ini(:)
	56	REAL , ALLOCATABLE :: dlon(:), dlat(:)
	57	REAL , ALLOCATABLE :: dlon_ini(:), dlat_ini(:)
	58	REAL , ALLOCATABLE :: champ_msk(:), champ(:)
	59	REAL , ALLOCATABLE :: work(:,:)
	60
	61	CHARACTER*25 title
	62
	63	C Declarations pour le champ interpole 2D
	64	REAL champint(iim,jjp1)
	65	real chmin,chmax
	66
	67	C Declarations pour le champ interpole 3D
	68	REAL champtime(iim,jjp1,tbid)
	69	REAL timeyear(tbid)
	70	REAL champan(iip1,jjp1,366)
	71
	72	C Declarations pour l'inteprolation verticale
	73	REAL ax(tbid), ay(tbid)
	74	REAL by
	75	REAL yder(tbid)
	76
	77
	78	INTEGER ierr
	79	INTEGER dimfirst(3)
	80	INTEGER dimlast(3)
	81	c
	82	INTEGER nid, ndim, ntim
	83	INTEGER dims(2), debut(2), epais(2)
	84	INTEGER id_tim
	85	INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
	86
	87	INTEGER i, j, k, l
	88	c Diverses variables locales
	89	REAL time
	90
	91	INTEGER longcles
	92	PARAMETER ( longcles = 20 )
	93	REAL clesphy0 ( longcles )
	94	#include "serre.h"
	95	INTEGER ncid,varid,ndimid(4),dimid
	96	character*30 namedim
	97
	98	pi = 4. * ATAN(1.)
	99	rad = 6 371 229.
	100	omeg = 4.* ASIN(1.)/(24.*3600.)
	101	g = 9.8
	102	daysec = 86400.
	103	kappa = 0.2857143
	104	cpp = 1004.70885
	105	dtvr = daysec/FLOAT(day_step)
	106	c
	107	C Traitement du relief au sol
	108	c
	109	write(,) 'Traitement du relief au sol pour fabriquer masque'
	110	ierr = NF_OPEN('Relief.nc', NF_NOWRITE, ncid)
	111
	112	if (ierr.ne.0) then
	113	print *, NF_STRERROR(ierr)
	114	STOP
	115	ENDIF
	116
	117	ierr = NF_INQ_VARID(ncid,'RELIEF',varid)
	118	if (ierr.ne.0) then
	119	print *, NF_STRERROR(ierr)
	120	STOP
	121	ENDIF
	122	ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
	123	if (ierr.ne.0) then
	124	print *, NF_STRERROR(ierr)
	125	STOP
	126	ENDIF
	127	ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
	128	if (ierr.ne.0) then
	129	print *, NF_STRERROR(ierr)
	130	STOP
	131	ENDIF
	132	print*,'variable ', namedim, 'dimension ', imdep
	133	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	134	if (ierr.ne.0) then
	135	print *, NF_STRERROR(ierr)
	136	STOP
	137	ENDIF
	138
	139	ALLOCATE( lonmsk_ini(imdep) )
	140	ALLOCATE( dlon_msk(imdep) )
	141
	142	#ifdef NC_DOUBLE
	143	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,lonmsk_ini)
	144	#else
	145	ierr = NF_GET_VAR_REAL(ncid,dimid,lonmsk_ini)
	146	#endif
	147
	148	c
	149	if (ierr.ne.0) then
	150	print *, NF_STRERROR(ierr)
	151	STOP
	152	ENDIF
	153	ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
	154	if (ierr.ne.0) then
	155	print *, NF_STRERROR(ierr)
	156	STOP
	157	ENDIF
	158	print*,'variable ', namedim, 'dimension ', jmdep
	159	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	160	if (ierr.ne.0) then
	161	print *, NF_STRERROR(ierr)
	162	STOP
	163	ENDIF
	164
	165	ALLOCATE( latmsk_ini(jmdep) )
	166	ALLOCATE( dlat_msk(jmdep) )
	167	ALLOCATE( champ_msk(imdep*jmdep) )
	168
	169	#ifdef NC_DOUBLE
	170	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,latmsk_ini)
	171	#else
	172	ierr = NF_GET_VAR_REAL(ncid,dimid,latmsk_ini)
	173	#endif
	174	c
	175	if (ierr.ne.0) then
	176	print *, NF_STRERROR(ierr)
	177	STOP
	178	ENDIF
	179	#ifdef NC_DOUBLE
	180	ierr = NF_GET_VAR_DOUBLE(ncid,varid,champ_msk)
	181	#else
	182	ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk)
	183	#endif
	184	c
	185	if (ierr.ne.0) then
	186	print *, NF_STRERROR(ierr)
	187	STOP
	188	ENDIF
	189	c
	190	title='RELIEF'
	191
	192	CALL conf_dat2d(title,imdep, jmdep, lonmsk_ini, latmsk_ini,
	193	. dlon_msk, dlat_msk, champ_msk, interbar )
	194
	195	CALL mask_c_o(imdep, jmdep, dlon_msk, dlat_msk,champ_msk,
	196	. iim, jjp1, rlonv, rlatu, champint)
	197	CALL gr_int_dyn(champint, masque, iim, jjp1)
	198	DO i = 1, iim
	199	masque(i,1) = FLOAT(NINT(masque(i,1)))
	200	masque(i,jjp1) = FLOAT(NINT(masque(i,jjp1)))
	201	ENDDO
	202	DO i = 1, iim
	203	DO j = 1, jjp1
	204	mask(i,j) = champint(i,j)
	205	ENDDO
	206	ENDDO
	207	CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0)
	208	ierr = NF_CLOSE(ncid)
	209	c
	210	c
	211	C Traitement de la rugosite
	212	c
	213	PRINT*, 'Traitement de la rugosite'
	214	ierr = NF_OPEN('Rugos.nc', NF_NOWRITE, ncid)
	215	if (ierr.ne.0) then
	216	print *, NF_STRERROR(ierr)
	217	STOP
	218	ENDIF
	219
	220	ierr = NF_INQ_VARID(ncid,'RUGOS',varid)
	221	if (ierr.ne.0) then
	222	print *, NF_STRERROR(ierr)
	223	STOP
	224	ENDIF
	225	ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
	226	if (ierr.ne.0) then
	227	print *, NF_STRERROR(ierr)
	228	STOP
	229	ENDIF
	230	ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
	231	if (ierr.ne.0) then
	232	print *, NF_STRERROR(ierr)
	233	STOP
	234	ENDIF
	235	print*,'variable ', namedim, 'dimension ', imdep
	236	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	237	if (ierr.ne.0) then
	238	print *, NF_STRERROR(ierr)
	239	STOP
	240	ENDIF
	241
	242	ALLOCATE( dlon_ini(imdep) )
	243	ALLOCATE( dlon(imdep) )
	244
	245	#ifdef NC_DOUBLE
	246	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
	247	#else
	248	ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
	249	#endif
	250	if (ierr.ne.0) then
	251	print *, NF_STRERROR(ierr)
	252	STOP
	253	ENDIF
	254	ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
	255	if (ierr.ne.0) then
	256	print *, NF_STRERROR(ierr)
	257	STOP
	258	ENDIF
	259	print*,'variable ', namedim, 'dimension ', jmdep
	260	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	261	if (ierr.ne.0) then
	262	print *, NF_STRERROR(ierr)
	263	STOP
	264	ENDIF
	265
	266	ALLOCATE( dlat_ini(jmdep) )
	267	ALLOCATE( dlat(jmdep) )
	268
	269	#ifdef NC_DOUBLE
	270	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
	271	#else
	272	ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
	273	#endif
	274	if (ierr.ne.0) then
	275	print *, NF_STRERROR(ierr)
	276	STOP
	277	ENDIF
	278	ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
	279	if (ierr.ne.0) then
	280	print *, NF_STRERROR(ierr)
	281	STOP
	282	ENDIF
	283	print*,'variable ', namedim, 'dimension ', lmdep
	284	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	285	if (ierr.ne.0) then
	286	print *, NF_STRERROR(ierr)
	287	STOP
	288	ENDIF
	289	#ifdef NC_DOUBLE
	290	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
	291	#else
	292	ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
	293	#endif
	294	if (ierr.ne.0) then
	295	print *, NF_STRERROR(ierr)
	296	STOP
	297	ENDIF
	298	c
	299	ALLOCATE( champ(imdep*jmdep) )
	300
	301	DO 200 l = 1, lmdep
	302	dimfirst(1) = 1
	303	dimfirst(2) = 1
	304	dimfirst(3) = l
	305	c
	306	dimlast(1) = imdep
	307	dimlast(2) = jmdep
	308	dimlast(3) = 1
	309	c
	310	PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
	311	print*,dimfirst,dimlast
	312	#ifdef NC_DOUBLE
	313	ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
	314	#else
	315	ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
	316	#endif
	317	if (ierr.ne.0) then
	318	print *, NF_STRERROR(ierr)
	319	STOP
	320	ENDIF
	321
	322	title = 'Rugosite Amip '
	323	c
	324	CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
	325	. dlon, dlat, champ, interbar )
	326
	327	IF ( interbar ) THEN
	328	DO j = 1, imdep * jmdep
	329	champ(j) = LOG(champ(j))
	330	ENDDO
	331
	332	IF( l.EQ.1 ) THEN
	333	WRITE(6,*) '-------------------------------------------------',
	334	,'------------------------'
	335	WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
	336	, ' pour la rugosite $$$ '
	337	WRITE(6,*) '-------------------------------------------------',
	338	,'------------------------'
	339	ENDIF
	340	CALL inter_barxy ( imdep,jmdep -1,dlon,dlat,champ ,
	341	, iim,jjm,rlonu,rlatv, jjp1,champint )
	342	DO j=1,jjp1
	343	DO i=1,iim
	344	champint(i,j)=EXP(champint(i,j))
	345	ENDDO
	346	ENDDO
	347
	348	DO j = 1, jjp1
	349	DO i = 1, iim
	350	IF(NINT(mask(i,j)).NE.1) THEN
	351	champint( i,j ) = 0.001
	352	ENDIF
	353	ENDDO
	354	ENDDO
	355	ELSE
	356	CALL rugosite(imdep, jmdep, dlon, dlat, champ,
	357	. iim, jjp1, rlonv, rlatu, champint, mask)
	358	ENDIF
	359	DO j = 1,jjp1
	360	DO i = 1, iim
	361	champtime (i,j,l) = champint(i,j)
	362	ENDDO
	363	ENDDO
	364	200 CONTINUE
	365	c
	366	DO l = 1, lmdep
	367	timeyear(l) = timecoord(l)
	368	ENDDO
	369
	370	PRINT 222, timeyear
	371	222 FORMAT(2x,' Time year ',10f6.1)
	372	c
	373
	374	PRINT*, 'Interpolation temporelle dans l annee'
	375
	376	DO j = 1, jjp1
	377	DO i = 1, iim
	378	DO l = 1, lmdep
	379	ax(l) = timeyear(l)
	380	ay(l) = champtime (i,j,l)
	381	ENDDO
	382	CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
	383	DO k = 1, 360
	384	time = FLOAT(k-1)
	385	CALL SPLINT(ax,ay,yder,lmdep,time,by)
	386	champan(i,j,k) = by
	387	ENDDO
	388	ENDDO
	389	ENDDO
	390	DO k = 1, 360
	391	DO j = 1, jjp1
	392	champan(iip1,j,k) = champan(1,j,k)
	393	ENDDO
	394	IF ( k.EQ.10 ) THEN
	395	DO j = 1, jjp1
	396	CALL minmax( iip1,champan(1,j,10),chmin,chmax )
	397	PRINT *,' Rugosite au temps 10 ', chmin,chmax,j
	398	ENDDO
	399	ENDIF
	400	ENDDO
	401	c
	402	DO k = 1, 360
	403	CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k), phy_rug(1,k))
	404	ENDDO
	405	c
	406	ierr = NF_CLOSE(ncid)
	407
	408	DEALLOCATE( dlon )
	409	DEALLOCATE( dlon_ini )
	410	DEALLOCATE( dlat )
	411	DEALLOCATE( dlat_ini )
	412	DEALLOCATE( champ )
	413	c
	414	c
	415	C Traitement de la glace oceanique
	416	c
	417	PRINT*, 'Traitement de la glace oceanique'
	418	ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
	419	if (ierr.ne.0) then
	420	print *, NF_STRERROR(ierr)
	421	STOP
	422	ENDIF
	423
	424	ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid)
	425	if (ierr.ne.0) then
	426	print *, NF_STRERROR(ierr)
	427	STOP
	428	ENDIF
	429	ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
	430	if (ierr.ne.0) then
	431	print *, NF_STRERROR(ierr)
	432	STOP
	433	ENDIF
	434	ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
	435	if (ierr.ne.0) then
	436	print *, NF_STRERROR(ierr)
	437	STOP
	438	ENDIF
	439	print*,'variable ', namedim, 'dimension ', imdep
	440	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	441	if (ierr.ne.0) then
	442	print *, NF_STRERROR(ierr)
	443	STOP
	444	ENDIF
	445
	446	ALLOCATE ( dlon_ini(imdep) )
	447	ALLOCATE ( dlon(imdep) )
	448
	449	#ifdef NC_DOUBLE
	450	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
	451	#else
	452	ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
	453	#endif
	454	if (ierr.ne.0) then
	455	print *, NF_STRERROR(ierr)
	456	STOP
	457	ENDIF
	458	ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
	459	if (ierr.ne.0) then
	460	print *, NF_STRERROR(ierr)
	461	STOP
	462	ENDIF
	463	print*,'variable ', namedim, jmdep
	464	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	465	if (ierr.ne.0) then
	466	print *, NF_STRERROR(ierr)
	467	STOP
	468	ENDIF
	469
	470	ALLOCATE ( dlat_ini(jmdep) )
	471	ALLOCATE ( dlat(jmdep) )
	472
	473	#ifdef NC_DOUBLE
	474	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
	475	#else
	476	ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
	477	#endif
	478	if (ierr.ne.0) then
	479	print *, NF_STRERROR(ierr)
	480	STOP
	481	ENDIF
	482	ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
	483	if (ierr.ne.0) then
	484	print *, NF_STRERROR(ierr)
	485	STOP
	486	ENDIF
	487	print*,'variable ', namedim, lmdep
	488	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	489	if (ierr.ne.0) then
	490	print *, NF_STRERROR(ierr)
	491	STOP
	492	ENDIF
	493	#ifdef NC_DOUBLE
	494	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
	495	#else
	496	ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
	497	#endif
	498	if (ierr.ne.0) then
	499	print *, NF_STRERROR(ierr)
	500	STOP
	501	ENDIF
	502	c
	503	ALLOCATE ( champ(imdep*jmdep) )
	504
	505	DO l = 1, lmdep
	506	dimfirst(1) = 1
	507	dimfirst(2) = 1
	508	dimfirst(3) = l
	509	c
	510	dimlast(1) = imdep
	511	dimlast(2) = jmdep
	512	dimlast(3) = 1
	513	c
	514	PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
	515	#ifdef NC_DOUBLE
	516	ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
	517	#else
	518	ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
	519	#endif
	520	if (ierr.ne.0) then
	521	print *, NF_STRERROR(ierr)
	522	STOP
	523	ENDIF
	524
	525	title = 'Sea-ice Amip '
	526	c
[425]	527	IF( oldice ) THEN
	528	CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
	529	. dlon, dlat, champ, .FALSE. )
	530	ELSE
	531	CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
[259]	532	. dlon, dlat, champ, interbar )
[425]	533	ENDIF
[259]	534	c
	535	IF( oldice ) THEN
	536	CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
	537	. iim, jjp1, rlonv, rlatu, champint )
	538	ELSEIF ( interbar ) THEN
	539	IF( l.EQ.1 ) THEN
	540	WRITE(6,*) '-------------------------------------------------',
	541	,'------------------------'
	542	WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
	543	, ' pour Sea-ice Amip $$$ '
	544	WRITE(6,*) '-------------------------------------------------',
	545	,'------------------------'
	546	ENDIF
	547
	548	CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
	549	, champ, iim, jjm, rlonu, rlatv, jjp1, champint )
	550	ELSE
	551	CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
	552	. iim, jjp1, rlonv, rlatu, champint )
	553	ENDIF
	554	DO j = 1,jjp1
	555	DO i = 1, iim
	556	champtime (i,j,l) = champint(i,j)
	557	ENDDO
	558	ENDDO
	559	ENDDO
	560	c
	561	DO l = 1, lmdep
	562	timeyear(l) = timecoord(l)
	563	ENDDO
	564	PRINT 222, timeyear
	565	c
	566	PRINT*, 'Interpolation temporelle'
	567	DO j = 1, jjp1
	568	DO i = 1, iim
	569	DO l = 1, lmdep
	570	ax(l) = timeyear(l)
	571	ay(l) = champtime (i,j,l)
	572	ENDDO
	573	CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
	574	DO k = 1, 360
	575	time = FLOAT(k-1)
	576	CALL SPLINT(ax,ay,yder,lmdep,time,by)
	577	champan(i,j,k) = by
	578	ENDDO
	579	ENDDO
	580	ENDDO
	581	DO k = 1, 360
	582	DO j = 1, jjp1
	583	champan(iip1, j, k) = champan(1, j, k)
	584	ENDDO
	585	IF ( k.EQ.10 ) THEN
	586	DO j = 1, jjp1
	587	CALL minmax( iip1,champan(1,j,10),chmin,chmax )
	588	PRINT *,' Sea ice au temps 10 ', chmin,chmax,j
	589	ENDDO
	590	ENDIF
	591	ENDDO
	592	c
	593
	594	DO k = 1, 360
	595	CALL gr_dyn_fi(1, iip1, jjp1, klon,
	596	. champan(1,1,k), phy_ice(1,k))
	597	DO i = 1, klon
	598	phy_nat(i,k) = phy_nat0(i)
	599	IF ( (phy_ice(i,k) - 0.5).GE.1.e-5 ) THEN
	600	IF (NINT(phy_nat0(i)).EQ.0) THEN
	601	phy_nat(i,k) = 3.0
	602	ELSE
	603	phy_nat(i,k) = 2.0
	604	ENDIF
	605	ENDIF
	606	IF( NINT(phy_nat(i,k)).EQ.0 ) THEN
	607	IF ( phy_rug(i,k).NE.0.001 ) phy_rug(i,k) = 0.001
	608	ENDIF
	609	ENDDO
	610
	611	ENDDO
	612	c
	613
	614	ierr = NF_CLOSE(ncid)
	615	c
	616	DEALLOCATE( dlon )
	617	DEALLOCATE( dlon_ini )
	618	DEALLOCATE( dlat )
	619	DEALLOCATE( dlat_ini )
	620	DEALLOCATE( champ )
	621
	622	477 continue
	623	c
	624	C Traitement de la sst
	625	c
	626	PRINT*, 'Traitement de la sst'
	627	ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
	628	if (ierr.ne.0) then
	629	print *, NF_STRERROR(ierr)
	630	STOP
	631	ENDIF
	632
	633	ierr = NF_INQ_VARID(ncid,'SST',varid)
	634	if (ierr.ne.0) then
	635	print *, NF_STRERROR(ierr)
	636	STOP
	637	ENDIF
	638	ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
	639	if (ierr.ne.0) then
	640	print *, NF_STRERROR(ierr)
	641	STOP
	642	ENDIF
	643	ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
	644	if (ierr.ne.0) then
	645	print *, NF_STRERROR(ierr)
	646	STOP
	647	ENDIF
	648	print*,'variable SST ', namedim,'dimension ', imdep
	649	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	650	if (ierr.ne.0) then
	651	print *, NF_STRERROR(ierr)
	652	STOP
	653	ENDIF
	654
	655	ALLOCATE( dlon_ini(imdep) )
	656	ALLOCATE( dlon(imdep) )
	657
	658	#ifdef NC_DOUBLE
	659	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
	660	#else
	661	ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
	662	#endif
	663
	664	if (ierr.ne.0) then
	665	print *, NF_STRERROR(ierr)
	666	STOP
	667	ENDIF
	668	ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
	669	if (ierr.ne.0) then
	670	print *, NF_STRERROR(ierr)
	671	STOP
	672	ENDIF
	673	print*,'variable SST ', namedim, 'dimension ', jmdep
	674	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	675	if (ierr.ne.0) then
	676	print *, NF_STRERROR(ierr)
	677	STOP
	678	ENDIF
	679
	680	ALLOCATE( dlat_ini(jmdep) )
	681	ALLOCATE( dlat(jmdep) )
	682
	683	#ifdef NC_DOUBLE
	684	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
	685	#else
	686	ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
	687	#endif
	688	if (ierr.ne.0) then
	689	print *, NF_STRERROR(ierr)
	690	STOP
	691	ENDIF
	692	ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
	693	if (ierr.ne.0) then
	694	print *, NF_STRERROR(ierr)
	695	STOP
	696	ENDIF
	697	print*,'variable ', namedim, 'dimension ', lmdep
	698	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	699	if (ierr.ne.0) then
	700	print *, NF_STRERROR(ierr)
	701	STOP
	702	ENDIF
	703	#ifdef NC_DOUBLE
	704	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
	705	#else
	706	ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
	707	#endif
	708	if (ierr.ne.0) then
	709	print *, NF_STRERROR(ierr)
	710	STOP
	711	ENDIF
	712
	713	ALLOCATE( champ(imdep*jmdep) )
	714	IF( extrap ) THEN
	715	ALLOCATE ( work(imdep,jmdep) )
	716	ENDIF
	717	c
	718	DO l = 1, lmdep
	719	dimfirst(1) = 1
	720	dimfirst(2) = 1
	721	dimfirst(3) = l
	722	c
	723	dimlast(1) = imdep
	724	dimlast(2) = jmdep
	725	dimlast(3) = 1
	726	c
	727	PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
	728	#ifdef NC_DOUBLE
	729	ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
	730	#else
	731	ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
	732	#endif
	733	if (ierr.ne.0) then
	734	print *, NF_STRERROR(ierr)
	735	STOP
	736	ENDIF
	737
	738	title='Sst Amip'
	739	c
	740	CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
	741	. dlon, dlat, champ, interbar )
	742	IF ( extrap ) THEN
	743	CALL extrapol(champ, imdep, jmdep, 999999.,.TRUE.,.TRUE.,2,work)
	744	ENDIF
	745	c
	746
	747	IF ( interbar ) THEN
	748	IF( l.EQ.1 ) THEN
	749	WRITE(6,*) '-------------------------------------------------',
	750	,'------------------------'
	751	WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
	752	, ' pour la Sst Amip $$$ '
	753	WRITE(6,*) '-------------------------------------------------',
	754	,'------------------------'
	755	ENDIF
	756	CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
	757	, champ, iim, jjm, rlonu, rlatv, jjp1, champint )
	758	ELSE
	759	CALL grille_m(imdep, jmdep, dlon, dlat, champ,
	760	. iim, jjp1, rlonv, rlatu, champint )
	761	ENDIF
	762
	763	DO j = 1,jjp1
	764	DO i = 1, iim
	765	champtime (i,j,l) = champint(i,j)
	766	ENDDO
	767	ENDDO
	768	ENDDO
	769	c
	770	DO l = 1, lmdep
	771	timeyear(l) = timecoord(l)
	772	ENDDO
	773	print 222, timeyear
	774	c
	775	C interpolation temporelle
	776	DO j = 1, jjp1
	777	DO i = 1, iim
	778	DO l = 1, lmdep
	779	ax(l) = timeyear(l)
	780	ay(l) = champtime (i,j,l)
	781	ENDDO
	782	CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
	783	DO k = 1, 360
	784	time = FLOAT(k-1)
	785	CALL SPLINT(ax,ay,yder,lmdep,time,by)
	786	champan(i,j,k) = by
	787	ENDDO
	788	ENDDO
	789	ENDDO
	790	DO k = 1, 360
	791	DO j = 1, jjp1
	792	champan(iip1,j,k) = champan(1,j,k)
	793	ENDDO
	794	IF ( k.EQ.10 ) THEN
	795	DO j = 1, jjp1
	796	CALL minmax( iip1,champan(1,j,10),chmin,chmax )
	797	PRINT *,' SST au temps 10 ', chmin,chmax,j
	798	ENDDO
	799	ENDIF
	800	ENDDO
	801	c
	802	DO k = 1, 360
	803	CALL gr_dyn_fi(1, iip1, jjp1, klon,
	804	. champan(1,1,k), phy_sst(1,k))
	805	ENDDO
	806	c
	807	ierr = NF_CLOSE(ncid)
	808	c
	809	DEALLOCATE( dlon )
	810	DEALLOCATE( dlon_ini )
	811	DEALLOCATE( dlat )
	812	DEALLOCATE( dlat_ini )
	813	DEALLOCATE( champ )
	814	c
	815	C Traitement de l'albedo
	816	c
	817	PRINT*, 'Traitement de l albedo'
	818	ierr = NF_OPEN('Albedo.nc', NF_NOWRITE, ncid)
	819	if (ierr.ne.0) then
	820	print *, NF_STRERROR(ierr)
	821	STOP
	822	ENDIF
	823	ierr = NF_INQ_VARID(ncid,'ALBEDO',varid)
	824	if (ierr.ne.0) then
	825	print *, NF_STRERROR(ierr)
	826	STOP
	827	ENDIF
	828	ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
	829	if (ierr.ne.0) then
	830	print *, NF_STRERROR(ierr)
	831	STOP
	832	ENDIF
	833	ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
	834	if (ierr.ne.0) then
	835	print *, NF_STRERROR(ierr)
	836	STOP
	837	ENDIF
	838	print*,'variable ', namedim, 'dimension ', imdep
	839	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	840	if (ierr.ne.0) then
	841	print *, NF_STRERROR(ierr)
	842	STOP
	843	ENDIF
	844
	845	ALLOCATE ( dlon_ini(imdep) )
	846	ALLOCATE ( dlon(imdep) )
	847
	848	#ifdef NC_DOUBLE
	849	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
	850	#else
	851	ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
	852	#endif
	853	if (ierr.ne.0) then
	854	print *, NF_STRERROR(ierr)
	855	STOP
	856	ENDIF
	857	ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
	858	if (ierr.ne.0) then
	859	print *, NF_STRERROR(ierr)
	860	STOP
	861	ENDIF
	862	print*,'variable ', namedim, 'dimension ', jmdep
	863	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	864	if (ierr.ne.0) then
	865	print *, NF_STRERROR(ierr)
	866	STOP
	867	ENDIF
	868
	869	ALLOCATE ( dlat_ini(jmdep) )
	870	ALLOCATE ( dlat(jmdep) )
	871
	872	#ifdef NC_DOUBLE
	873	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
	874	#else
	875	ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
	876	#endif
	877	if (ierr.ne.0) then
	878	print *, NF_STRERROR(ierr)
	879	STOP
	880	ENDIF
	881	ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
	882	if (ierr.ne.0) then
	883	print *, NF_STRERROR(ierr)
	884	STOP
	885	ENDIF
	886	print*,'variable ', namedim, 'dimension ', lmdep
	887	ierr = NF_INQ_VARID(ncid,namedim,dimid)
	888	if (ierr.ne.0) then
	889	print *, NF_STRERROR(ierr)
	890	STOP
	891	ENDIF
	892	#ifdef NC_DOUBLE
	893	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
	894	#else
	895	ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
	896	#endif
	897	if (ierr.ne.0) then
	898	print *, NF_STRERROR(ierr)
	899	STOP
	900	ENDIF
	901	c
	902	ALLOCATE ( champ(imdep*jmdep) )
	903
	904	DO l = 1, lmdep
	905	dimfirst(1) = 1
	906	dimfirst(2) = 1
	907	dimfirst(3) = l
	908	c
	909	dimlast(1) = imdep
	910	dimlast(2) = jmdep
	911	dimlast(3) = 1
	912	c
	913	PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
	914	#ifdef NC_DOUBLE
	915	ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
	916	#else
	917	ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
	918	#endif
	919	if (ierr.ne.0) then
	920	print *, NF_STRERROR(ierr)
	921	STOP
	922	ENDIF
	923
	924	title='Albedo Amip'
	925	c
	926	CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
	927	. dlon, dlat, champ, interbar )
	928	c
	929	c
	930	IF ( interbar ) THEN
	931	IF( l.EQ.1 ) THEN
	932	WRITE(6,*) '-------------------------------------------------',
	933	,'------------------------'
	934	WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
	935	, ' pour l Albedo Amip $$$ '
	936	WRITE(6,*) '-------------------------------------------------',
	937	,'------------------------'
	938	ENDIF
	939
	940	CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
	941	, champ, iim, jjm, rlonu, rlatv, jjp1, champint )
	942	ELSE
	943	CALL grille_m(imdep, jmdep, dlon, dlat, champ,
	944	. iim, jjp1, rlonv, rlatu, champint )
	945	ENDIF
	946	c
	947	DO j = 1,jjp1
	948	DO i = 1, iim
	949	champtime (i, j, l) = champint(i, j)
	950	ENDDO
	951	ENDDO
	952	ENDDO
	953	c
	954	DO l = 1, lmdep
	955	timeyear(l) = timecoord(l)
	956	ENDDO
	957	print 222, timeyear
	958	c
	959	C interpolation temporelle
	960	DO j = 1, jjp1
	961	DO i = 1, iim
	962	DO l = 1, lmdep
	963	ax(l) = timeyear(l)
	964	ay(l) = champtime (i, j, l)
	965	ENDDO
	966	CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
	967	DO k = 1, 360
	968	time = FLOAT(k-1)
	969	CALL SPLINT(ax,ay,yder,lmdep,time,by)
	970	champan(i,j,k) = by
	971	ENDDO
	972	ENDDO
	973	ENDDO
	974	DO k = 1, 360
	975	DO j = 1, jjp1
	976	champan(iip1, j, k) = champan(1, j, k)
	977	ENDDO
	978	IF ( k.EQ.10 ) THEN
	979	DO j = 1, jjp1
	980	CALL minmax( iip1,champan(1,j,10),chmin,chmax )
	981	PRINT *,' Albedo au temps 10 ', chmin,chmax,j
	982	ENDDO
	983	ENDIF
	984	ENDDO
	985	c
	986	DO k = 1, 360
	987	CALL gr_dyn_fi(1, iip1, jjp1, klon,
	988	. champan(1,1,k), phy_alb(1,k))
	989	ENDDO
	990	c
	991	ierr = NF_CLOSE(ncid)
	992	c
	993	c
	994	DO k = 1, 360
	995	DO i = 1, klon
	996	phy_bil(i,k) = 0.0
	997	ENDDO
	998	ENDDO
	999	c
	1000	PRINT*, 'Ecriture du fichier limit'
	1001	c
	1002	ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
	1003	c
	1004	ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
	1005	. "Fichier conditions aux limites")
	1006	ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
	1007	ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
	1008	c
	1009	dims(1) = ndim
	1010	dims(2) = ntim
	1011	c
	1012	ccc ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
	1013	ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
	1014	ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
	1015	. "Jour dans l annee")
	1016	ccc ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
	1017	ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
	1018	ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
	1019	. "Nature du sol (0,1,2,3)")
	1020	ccc ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
	1021	ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
	1022	ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
	1023	. "Temperature superficielle de la mer")
	1024	ccc ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
	1025	ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
	1026	ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
	1027	. "Reference flux de chaleur au sol")
	1028	ccc ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
	1029	ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
	1030	ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
	1031	. "Albedo a la surface")
	1032	ccc ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
	1033	ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
	1034	ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
	1035	. "Rugosite")
	1036	c
	1037	ierr = NF_ENDDEF(nid)
	1038	c
	1039	DO k = 1, 360
	1040	c
	1041	debut(1) = 1
	1042	debut(2) = k
	1043	epais(1) = klon
	1044	epais(2) = 1
	1045	c
	1046	#ifdef NC_DOUBLE
	1047	ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
	1048	ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
	1049	ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
	1050	ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
	1051	ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
	1052	ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
	1053	#else
	1054	ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
	1055	ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
	1056	ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
	1057	ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
	1058	ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
	1059	ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
	1060	#endif
	1061	c
	1062	ENDDO
	1063	c
	1064	ierr = NF_CLOSE(nid)
	1065	c
	1066	STOP
	1067	END

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

Download in other formats:

Original Format