Context Navigation

create_limit.F @ 278

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

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

Download in other formats:

Original Format