Context Navigation

etat0_netcdf.F @ 359

Last change on this file since 359 was 359, checked in by lmdzadmin, 22 years ago
Menage sur les etats initiaux LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 21.4 KB

Line
1	c
2	c $Header$
3	c
4	SUBROUTINE etat0_netcdf (interbar, masque, pctsrf)
5
6	USE startvar
7	USE ioipsl
8	!
9	IMPLICIT NONE
10	!
11	#include "netcdf.inc"
12	#include "dimensions.h"
13	#include "paramet.h"
14	!
15	!
16	! INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2,
17	! .KLON=KFDIA-KIDIA+1,KLEV=llm
18	!
19	#include "comgeom2.h"
20	#include "comvert.h"
21	#include "comconst.h"
22	#include "indicesol.h"
23	#include "dimphy.h"
24	#include "dimsoil.h"
25	#include "temps.h"
26	!
27	LOGICAL interbar
28	REAL :: latfi(klon), lonfi(klon)
29	REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
30	. psol(iip1, jjp1), phis(iip1, jjp1)
31	REAL :: p3d(iip1, jjp1, llm+1)
32	REAL :: uvent(iip1, jjp1, llm)
33	REAL :: vvent(iip1, jjm, llm)
34	REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
35	REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
36	REAL :: tsol(klon), qsol(klon), sn(klon)
37	REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf)
38	REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
39	REAL :: tsoil(klon,nsoilmx,nbsrf)
40	REAL :: radsol(klon),rain_fall(klon), snow_fall(klon)
41	REAL :: solsw(klon), sollw(klon), fder(klon)
42	REAL :: deltat(klon), frugs(klon,nbsrf), agesno(klon,nbsrf)
43	REAL :: rugmer(klon)
44	REAL :: zmea(iip1jjp1), zstd(iip1jjp1)
45	REAL :: zsig(iip1jjp1), zgam(iip1jjp1), zthe(iip1*jjp1)
46	REAL :: zpic(iip1jjp1), zval(iip1jjp1), rugsrel(iip1*jjp1)
47	REAL :: qd(iip1, jjp1, llm)
48	REAL :: pctsrf(klon, nbsrf)
49	REAL :: t_ancien(klon,klev), q_ancien(klon,klev) !
50	! declarations pour lecture glace de mer
51	REAL :: rugv(klon)
52	INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
53	INTEGER :: itaul(1), fid
54	REAL :: lev(1), date
55	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
56	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic
57	REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
58	REAL :: flic_tmp(iip1, jjp1)
59	REAL :: champint(iim, jjp1)
60	!
61
62	CHARACTER*80 :: varname
63	!
64	INTEGER :: i,j, ig, l, ji,ii1,ii2
65	REAL :: xpi
66	!
67	REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
68	REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
69	REAL :: workvar(iip1,jjp1,llm)
70	!
71	REAL :: prefkap, unskap
72	!
73	REAL :: q_sat
74	EXTERNAL q_sat
75	real :: time_step,t_ops,t_wrt
76
77	#include "comdissnew.h"
78	#include "control.h"
79	#include "serre.h"
80	#include "clesphys.h"
81
82	INTEGER :: longcles
83	PARAMETER ( longcles = 20 )
84	REAL :: clesphy0 ( longcles )
85	REAL :: p(iip1,jjp1,llm)
86	INTEGER :: itau, iday
87	REAL :: masse(iip1,jjp1,llm)
88	REAL :: xpn,xps,xppn(iim),xpps(iim)
89	real :: time
90	REAL :: phi(ip1jmp1,llm)
91	REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
92	REAL :: w(ip1jmp1,llm)
93	REAL ::phystep,co2_ppm,solaire
94	INTEGER :: radpas
95	real zrel(iip1*jjp1),chmin,chmax
96
97	CHARACTER*80 :: visu_file
98	INTEGER :: visuid
99
100	! pour la lecture du fichier masque ocean
101	integer :: nid_o2a
102	logical :: couple = .false.
103	INTEGER :: iml_omask, jml_omask
104	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
105	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask
106	REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
107	real, dimension(klon) :: ocemask_fi
108	integer :: isst(klon-2)
109	real zx_tmp_2d(iim,jjp1)
110	!
111	! Constantes
112	!
113	pi = 4. * ATAN(1.)
114	rad = 6371229.
115	omeg = 4.* ASIN(1.)/(24.*3600.)
116	g = 9.8
117	daysec = 86400.
118	kappa = 0.2857143
119	cpp = 1004.70885
120	!
121	preff = 101325.
122	unskap = 1./kappa
123	!
124	jmp1 = jjm + 1
125	!
126	! Construct a grid
127	!
128
129	! CALL defrun_new(99,.TRUE.,clesphy0)
130	CALL conf_gcm( 99, .TRUE. , clesphy0 )
131
132	dtvr = daysec/FLOAT(day_step)
133	print*,'dtvr',dtvr
134
135	CALL inicons0()
136	CALL inigeom()
137	!
138	CALL inifilr()
139	!
140	latfi(1) = ASIN(1.0)
141	DO j = 2, jjm
142	DO i = 1, iim
143	latfi((j-2)*iim+1+i)= rlatu(j)
144	ENDDO
145	ENDDO
146	latfi(klon) = - ASIN(1.0)
147	!
148	lonfi(1) = 0.0
149	DO j = 2, jjm
150	DO i = 1, iim
151	lonfi((j-2)*iim+1+i) = rlonv(i)
152	ENDDO
153	ENDDO
154	lonfi(klon) = 0.0
155	!
156	xpi = 2.0 * ASIN(1.0)
157	DO ig = 1, klon
158	latfi(ig) = latfi(ig) * 180.0 / xpi
159	lonfi(ig) = lonfi(ig) * 180.0 / xpi
160	ENDDO
161	!
162	varname = 'relief'
163	! This line needs to be replaced by a call to restget to get the values in the restart file
164	orog(:,:) = 0.0
165	CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
166	, jjm ,rlonu,rlatv , interbar )
167	!
168	WRITE(,) 'OUT OF GET VARIABLE : Relief'
169	! WRITE(*,'(49I1)') INT(orog(:,:))
170	!
171	varname = 'rugosite'
172	! This line needs to be replaced by a call to restget to get the values in the restart file
173	rugo(:,:) = 0.0
174	CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
175	, jjm, rlonu,rlatv , interbar )
176	!
177	WRITE(,) 'OUT OF GET VARIABLE : Rugosite'
178	! WRITE(,'(49I1)') INT(rugo(:,:)10)
179	!
180	varname = 'masque'
181	! This line needs to be replaced by a call to restget to get the values in the restart file
182	masque(:,:) = 0.0
183
184	CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
185	, jjm ,rlonu,rlatv , interbar )
186	! masque = nint(masque)
187	!
188	WRITE(,) 'MASQUE construit : Masque'
189	WRITE(*,'(97I1)') nINT(masque(:,:))
190	!
191	!
192
193
194	C
195	C on initialise les sous surfaces
196	C
197	pctsrf=0.
198	!cree le masque a partir du fichier relief
199	call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
200	WHERE (zmasq(1 : klon) .LT. EPSFRA)
201	zmasq(1 : klon) = 0.
202	END WHERE
203	WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
204	zmasq(1 : klon) = 1.
205	END WHERE
206	c
207	varname = 'psol'
208	psol(:,:) = 0.0
209	CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
210	, jjm ,rlonu,rlatv , interbar )
211	!
212	! Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
213	! anyway.
214	!
215	! WRITE(,) 'PSOL :', psol(10,20)
216	! WRITE(,) ap(:), bp(:)
217	CALL pression(ip1jmp1, ap, bp, psol, p3d)
218	! WRITE(,) 'P3D :', p3d(10,20,:)
219	CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
220	! WRITE(,) 'PK:', pk(10,20,:)
221	!
222	!
223	!
224	prefkap = preff ** kappa
225	! WRITE(,) 'unskap, cpp, preff :', unskap, cpp, preff
226	DO l = 1, llm
227	DO j=1,jjp1
228	DO i =1, iip1
229	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
230	ENDDO
231	ENDDO
232	ENDDO
233	!
234	! WRITE(,) 'PLS :', pls(10,20,:)
235	!
236	varname = 'surfgeo'
237	phis(:,:) = 0.0
238	CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
239	, jjm ,rlonu,rlatv, interbar )
240	!
241	varname = 'u'
242	uvent(:,:,:) = 0.0
243	CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
244	. workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
245	!
246	varname = 'v'
247	vvent(:,:,:) = 0.0
248	CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
249	. workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
250	!
251	varname = 't'
252	t3d(:,:,:) = 0.0
253	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
254	. workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
255	!
256	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
257	. maxval(t3d(:,:,:))
258	varname = 'tpot'
259	tpot(:,:,:) = 0.0
260	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
261	. pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
262	!
263	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
264	. maxval(t3d(:,:,:))
265	WRITE(,) 'PLS min,max:',minval(pls(:,:,:)),
266	. maxval(pls(:,:,:))
267	DO l = 1, llm
268	DO j=1,jjp1
269	DO i =1, iip1-1
270	qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
271	ENDDO
272	qsat(iip1,j,l) = qsat(1,j,l)
273	ENDDO
274	ENDDO
275	WRITE(,) 'QSAT min,max:',minval(qsat(:,:,:)),
276	. maxval(qsat(:,:,:))
277	!
278	WRITE(,) 'QSAT :', qsat(10,20,:)
279	!
280	varname = 'q'
281	qd(:,:,:) = 0.0
282	q3d(:,:,:,:) = 0.0
283	WRITE(,) 'QSAT min,max:',minval(qsat(:,:,:)),
284	. maxval(qsat(:,:,:))
285	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
286	. qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
287	q3d(:,:,:,1) = qd(:,:,:)
288	!
289	varname = 'tsol'
290	! This line needs to be replaced by a call to restget to get the values in the restart file
291	tsol(:) = 0.0
292	CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
293	. jjm, rlonu, rlatv , interbar )
294	!
295	WRITE(,) 'TSOL construit :'
296	! WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
297	!
298	varname = 'qsol'
299	qsol(:) = 0.0
300	CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0,
301	. jjm, rlonu, rlatv , interbar )
302	!
303	varname = 'snow'
304	sn(:) = 0.0
305	CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0,
306	. jjm, rlonu, rlatv , interbar )
307	!
308	varname = 'rads'
309	radsol(:) = 0.0
310	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
311	. jjm, rlonu, rlatv , interbar )
312	!
313	varname = 'deltat'
314	deltat(:) = 0.0
315	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0,
316	. jjm, rlonu, rlatv , interbar )
317	!
318	varname = 'rugmer'
319	rugmer(:) = 0.0
320	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
321	. jjm, rlonu, rlatv , interbar )
322	!
323	! varname = 'agesno'
324	! agesno(:) = 0.0
325	! CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
326	! . jjm, rlonu, rlatv , interbar )
327
328	varname = 'zmea'
329	zmea(:) = 0.0
330	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
331	. jjm, rlonu, rlatv , interbar )
332
333	varname = 'zstd'
334	zstd(:) = 0.0
335	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
336	. jjm, rlonu, rlatv , interbar )
337	varname = 'zsig'
338	zsig(:) = 0.0
339	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
340	. jjm, rlonu, rlatv , interbar )
341	varname = 'zgam'
342	zgam(:) = 0.0
343	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
344	. jjm, rlonu, rlatv , interbar )
345	varname = 'zthe'
346	zthe(:) = 0.0
347	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
348	. jjm, rlonu, rlatv , interbar )
349	varname = 'zpic'
350	zpic(:) = 0.0
351	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
352	. jjm, rlonu, rlatv , interbar )
353	varname = 'zval'
354	zval(:) = 0.0
355	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
356	. jjm, rlonu, rlatv , interbar )
357	c
358	rugsrel(:) = 0.0
359	IF(ok_orodr) THEN
360	DO i = 1, iip1* jjp1
361	rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
362	ENDDO
363	ENDIF
364	C
365	C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
366	C utilise pour calculer les poids et pour assurer l'adequation entre les
367	C fractions d'ocean vu par l'atmosphere et l'ocean
368	C
369
370	write(,)'Essai de lecture masque ocean'
371	iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
372	if (iret .ne. 0) then
373	write(,)'ATTENTION!! pas de fichier o2a.nc trouve'
374	write(,)'Run force'
375	else
376	couple = .true.
377	iret = nf_close(nid_o2a)
378	call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
379	$ , nid_o2a)
380	if (iml_omask /= iim .or. jml_omask /= jjp1) then
381	write(,)'Dimensions non compatibles pour masque ocean'
382	write(,)'iim = ',iim,' iml_omask = ',iml_omask
383	write(,)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
384	stop
385	endif
386	ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
387	ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
388	ALLOCATE(dlon_omask(iml_omask), stat=iret)
389	ALLOCATE(dlat_omask(jml_omask), stat=iret)
390	ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
391	ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
392	CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
393	$ , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
394	CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
395	$ ttm_tmp, 1, 1, ocetmp)
396	CALL flinclo(fid)
397	dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
398	dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
399	ocemask = ocetmp
400	if (dlat_omask(1) < dlat_omask(jml_omask)) then
401	do j = 1, jml_omask
402	ocemask(:,j) = ocetmp(:,jml_omask-j+1)
403	enddo
404	endif
405	C
406	C passage masque ocean a la grille physique
407	C
408
409	write(,)'ocemask '
410	write(*,'(96i1)')int(ocemask)
411	ocemask_fi(1) = ocemask(1,1)
412	do j = 2, jjm
413	do i = 1, iim
414	ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
415	enddo
416	enddo
417	ocemask_fi(klon) = ocemask(1,jjp1)
418	zmasq = 1. - ocemask_fi
419	isst = 0
420	where (ocemask_fi(2:klon-1) >0.) isst = 1
421	endif
422
423
424	C
425	C lecture du fichier glace de terre pour fixer la fraction de terre
426	C et de glace de terre
427	C
428	CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
429	$ , fid)
430	ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
431	ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
432	ALLOCATE(dlon_lic(iml_lic), stat=iret)
433	ALLOCATE(dlat_lic(jml_lic), stat=iret)
434	ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
435	CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
436	$ , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
437	CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
438	$ , 1, 1, fraclic)
439	CALL flinclo(fid)
440	C
441	C interpolation sur la grille T du modele
442	C
443	WRITE(,) 'dimensions de landice iml_lic, jml_lic : ',
444	$ iml_lic, jml_lic
445	c
446	C sil les coordonnees sont en degres, on les transforme
447	C
448	IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) ) THEN
449	lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
450	ENDIF
451	IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
452	lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
453	ENDIF
454
455	dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
456	dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
457	C
458	CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
459	$ ,iim, jjp1,
460	$ rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
461	c$$$ flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
462	flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
463	C
464	C passage sur la grille physique
465	C
466	CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
467	$ pctsrf(1:klon, is_lic))
468	C adequation avec le maque terre/mer
469	c zmasq(157) = 0.
470	WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
471	pctsrf(1 : klon, is_lic) = 0.
472	END WHERE
473	WHERE (zmasq( 1 : klon) .LT. EPSFRA)
474	pctsrf(1 : klon, is_lic) = 0.
475	END WHERE
476	pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
477	DO ji = 1, klon
478	IF (zmasq(ji) .GT. EPSFRA) THEN
479	IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
480	pctsrf(ji, is_lic) = zmasq(ji)
481	pctsrf(ji, is_ter) = 0.
482	ELSE
483	pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
484	IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
485	pctsrf(ji,is_ter) = 0.
486	pctsrf(ji, is_lic) = zmasq(ji)
487	ENDIF
488	ENDIF
489	ENDIF
490	END DO
491	C
492	C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
493	C
494	pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
495
496
497	WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
498	pctsrf(1 : klon, is_oce) = 0.
499	END WHERE
500
501	if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
502
503	isst = 0
504	where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
505	C
506	C verif que somme des sous surface = 1
507	C
508	ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
509	$ .GT. EPSFRA)
510	IF (ji .NE. 0) THEN
511	WRITE(,) 'pb repartition sous maille pour ',ji,' points'
512	ENDIF
513
514	! where (pctsrf(1:klon, is_ter) >= .5)
515	! pctsrf(1:klon, is_ter) = 1.
516	! pctsrf(1:klon, is_oce) = 0.
517	! pctsrf(1:klon, is_sic) = 0.
518	! pctsrf(1:klon, is_lic) = 0.
519	! zmasq = 1.
520	! endwhere
521	! where (pctsrf(1:klon, is_lic) >= .5)
522	! pctsrf(1:klon, is_ter) = 0.
523	! pctsrf(1:klon, is_oce) = 0.
524	! pctsrf(1:klon, is_sic) = 0.
525	! pctsrf(1:klon, is_lic) = 1.
526	! zmasq = 1.
527	! endwhere
528	! where (pctsrf(1:klon, is_oce) >= .5)
529	! pctsrf(1:klon, is_ter) = 0.
530	! pctsrf(1:klon, is_oce) = 1.
531	! pctsrf(1:klon, is_sic) = 0.
532	! pctsrf(1:klon, is_lic) = 0.
533	! zmasq = 0.
534	! endwhere
535	! where (pctsrf(1:klon, is_sic) >= .5)
536	! pctsrf(1:klon, is_ter) = 0.
537	! pctsrf(1:klon, is_oce) = 0.
538	! pctsrf(1:klon, is_sic) = 1.
539	! pctsrf(1:klon, is_lic) = 0.
540	! zmasq = 0.
541	! endwhere
542	! call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
543	C
544	C verif que somme des sous surface = 1
545	C
546	! ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
547	! $ .GT. EPSFRA)
548	! IF (ji .NE. 0) THEN
549	! WRITE(,) 'pb repartition sous maille pour ',ji,' points'
550	! ENDIF
551
552	CALL gr_fi_ecrit(1,klon,iim,jjp1,zmasq,zx_tmp_2d)
553	write(,)'zmasq = '
554	write(*,'(96i1)')nint(zx_tmp_2d)
555	call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
556	WRITE(,) 'MASQUE construit : Masque'
557	WRITE(*,'(97I1)') nINT(masque(:,:))
558
559
560
561	C Calcul intermediaire
562	c
563	CALL massdair( p3d, masse )
564	c
565
566	print *,' ALPHAX ',alphax
567
568	DO l = 1, llm
569	DO i = 1, iim
570	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
571	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
572	ENDDO
573	xpn = SUM(xppn)/apoln
574	xps = SUM(xpps)/apols
575	DO i = 1, iip1
576	masse( i , 1 , l ) = xpn
577	masse( i , jjp1 , l ) = xps
578	ENDDO
579	ENDDO
580	q3d(iip1,:,:,:) = q3d(1,:,:,:)
581	phis(iip1,:) = phis(1,:)
582
583	C Ecriture
584
585
586	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
587	* tetagdiv, tetagrot , tetatemp )
588	print*,'sortie inidissip'
589	itau = 0
590	itau_dyn = 0
591	itau_phy = 0
592	iday = dayref +itau/day_step
593	time = FLOAT(itau-(iday-dayref)*day_step)/day_step
594	c
595	IF(time.GT.1) THEN
596	time = time - 1
597	iday = iday + 1
598	ENDIF
599	day_ref = dayref
600	annee_ref = anneeref
601
602	CALL geopot ( ip1jmp1, tpot , pk , pks, phis , phi )
603	print*,'sortie geopot'
604
605	CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
606	* phi,w, pbaru,pbarv,time+iday-dayref )
607	print*,'sortie caldyn0'
608	CALL dynredem0("start.nc",dayref,phis,nqmx)
609	print*,'sortie dynredem0'
610	CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
611	. psol)
612	print*,'sortie dynredem1'
613	C
614	C Ecriture etat initial physique
615	C
616	write(,)'phystep ',dtvr,iphysiq,nbapp_rad
617	phystep = dtvr * FLOAT(iphysiq)
618	radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) )
619	write(,)'phystep =', phystep, radpas
620	co2_ppm = 330.0
621	solaire = 1370.0
622
623	c
624	c Initialisation
625	c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
626	c
627	tsolsrf(:,is_ter) = tsol
628	tsolsrf(:,is_lic) = tsol
629	tsolsrf(:,is_oce) = tsol
630	tsolsrf(:,is_sic) = tsol
631	snsrf(:,is_ter) = sn
632	snsrf(:,is_lic) = sn
633	snsrf(:,is_oce) = sn
634	snsrf(:,is_sic) = sn
635	albe(:,is_ter) = 0.08
636	albe(:,is_lic) = 0.6
637	albe(:,is_oce) = 0.5
638	albe(:,is_sic) = 0.6
639	evap(:,:) = 0.
640	qsolsrf(:,is_ter) = 150
641	qsolsrf(:,is_lic) = 150
642	qsolsrf(:,is_oce) = 150.
643	qsolsrf(:,is_sic) = 150.
644	do i = 1, nbsrf
645	do j = 1, nsoilmx
646	tsoil(:,j,i) = tsol
647	enddo
648	enddo
649	rain_fall = 0.; snow_fall = 0.
650	solsw = 165.
651	sollw = -53.
652	t_ancien = 273.15
653	q_ancien = 0.
654	agesno = 0.
655	deltat = 0.
656	frugs(1:klon,is_oce) = rugmer(1:klon)
657	frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
658	frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
659	frugs(1:klon,is_sic) = 0.001
660	fder = 0.0
661
662	call phyredem("startphy.nc",phystep,radpas, co2_ppm, solaire,
663	$ latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, snsrf,
664	$ albe, evap, rain_fall, snow_fall, solsw, sollw, fder,
665	$ radsol, frugs, agesno,
666	$ zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
667	$ t_ancien, q_ancien)
668
669	C Sortie Visu pour les champs dynamiques
670	print*,'sortie visu'
671	time_step = 1.
672	t_ops = 2.
673	t_wrt = 2.
674	itau = 2.
675	visu_file='Etat0_visu.nc'
676	CALL initdynav(visu_file,dayref,anneeref,time_step,
677	. t_ops, t_wrt, nqmx, visuid)
678	CALL writedynav(visuid, nqmx, itau,vvent ,
679	. uvent,tpot,pk,phi,q3d,masse,psol,phis)
680	print*,'entree histclo'
681	CALL histclo
682	RETURN
683	!
684	END SUBROUTINE etat0_netcdf

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

Download in other formats:

Original Format