Context Navigation

etat0_netcdf.F @ 1270

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

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

Download in other formats:

Original Format