Context Navigation

create_limit.F @ 495

Last change on this file since 495 was 326, checked in by lmdzadmin, 23 years ago
Pour etre compatible avec le startget de LeVan? LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 32.7 KB

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

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

Download in other formats:

Original Format