Context Navigation

etat0_netcdf.F @ 802

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

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

Download in other formats:

Original Format