Context Navigation

create_limit.F @ 228

Last change on this file since 228 was 177, checked in by lmdzadmin, 23 years ago

Lots of stuff, plus particulierement:

appel a ORCHIDEE en etat de marche (pb de grille subsiste)
modifs de Pascale sur soil dans le cas ou ok_veget=false -

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

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

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

Download in other formats:

Original Format