Context Navigation

etat0_netcdf.F @ 1189

Last change on this file since 1189 was 630, checked in by Laurent Fairhead, 20 years ago
Import d'une version parallele de la dynamique YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 21.8 KB

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

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

Download in other formats:

Original Format