Context Navigation

etat0_netcdf.F @ 354

Last change on this file since 354 was 354, checked in by lmdzadmin, 23 years ago
include du mauvais fichier apres un peu de menage LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 21.0 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
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	!
187	WRITE(,) 'MASQUE construit : Masque'
188	! WRITE(*,'(97I1)') nINT(masque(:,:))
189	!
190	!
191
192
193	C
194	C on initialise les sous surfaces
195	C
196	pctsrf=0.
197	!cree le masque a partir du fichier relief
198	call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
199	WHERE (zmasq(1 : klon) .LT. EPSFRA)
200	zmasq(1 : klon) = 0.
201	END WHERE
202	WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
203	zmasq(1 : klon) = 1.
204	END WHERE
205
206
207
208
209	varname = 'psol'
210	psol(:,:) = 0.0
211	CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
212	, jjm ,rlonu,rlatv , interbar )
213	!
214	! Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
215	! anyway.
216	!
217	! WRITE(,) 'PSOL :', psol(10,20)
218	! WRITE(,) ap(:), bp(:)
219	CALL pression(ip1jmp1, ap, bp, psol, p3d)
220	! WRITE(,) 'P3D :', p3d(10,20,:)
221	CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
222	! WRITE(,) 'PK:', pk(10,20,:)
223	!
224	!
225	!
226	prefkap = preff ** kappa
227	! WRITE(,) 'unskap, cpp, preff :', unskap, cpp, preff
228	DO l = 1, llm
229	DO j=1,jjp1
230	DO i =1, iip1
231	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
232	ENDDO
233	ENDDO
234	ENDDO
235	!
236	! WRITE(,) 'PLS :', pls(10,20,:)
237	!
238	varname = 'surfgeo'
239	phis(:,:) = 0.0
240	CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
241	, jjm ,rlonu,rlatv, interbar )
242	!
243	varname = 'u'
244	uvent(:,:,:) = 0.0
245	CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
246	. workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
247	!
248	varname = 'v'
249	vvent(:,:,:) = 0.0
250	CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
251	. workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
252	!
253	varname = 't'
254	t3d(:,:,:) = 0.0
255	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
256	. workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
257	!
258	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
259	. maxval(t3d(:,:,:))
260	varname = 'tpot'
261	tpot(:,:,:) = 0.0
262	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
263	. pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
264	!
265	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
266	. maxval(t3d(:,:,:))
267	WRITE(,) 'PLS min,max:',minval(pls(:,:,:)),
268	. maxval(pls(:,:,:))
269	DO l = 1, llm
270	DO j=1,jjp1
271	DO i =1, iip1-1
272	qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
273	ENDDO
274	qsat(iip1,j,l) = qsat(1,j,l)
275	ENDDO
276	ENDDO
277	WRITE(,) 'QSAT min,max:',minval(qsat(:,:,:)),
278	. maxval(qsat(:,:,:))
279	!
280	WRITE(,) 'QSAT :', qsat(10,20,:)
281	!
282	varname = 'q'
283	qd(:,:,:) = 0.0
284	q3d(:,:,:,:) = 0.0
285	WRITE(,) 'QSAT min,max:',minval(qsat(:,:,:)),
286	. maxval(qsat(:,:,:))
287	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
288	. qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
289	q3d(:,:,:,1) = qd(:,:,:)
290	!
291	varname = 'tsol'
292	! This line needs to be replaced by a call to restget to get the values in the restart file
293	tsol(:) = 0.0
294	CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
295	. jjm, rlonu, rlatv , interbar )
296	!
297	WRITE(,) 'TSOL construit :'
298	! WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
299	!
300	varname = 'qsol'
301	qsol(:) = 0.0
302	CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0,
303	. jjm, rlonu, rlatv , interbar )
304	!
305	varname = 'snow'
306	sn(:) = 0.0
307	CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0,
308	. jjm, rlonu, rlatv , interbar )
309	!
310	varname = 'rads'
311	radsol(:) = 0.0
312	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
313	. jjm, rlonu, rlatv , interbar )
314	!
315	varname = 'deltat'
316	deltat(:) = 0.0
317	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0,
318	. jjm, rlonu, rlatv , interbar )
319	!
320	varname = 'rugmer'
321	rugmer(:) = 0.0
322	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
323	. jjm, rlonu, rlatv , interbar )
324	!
325	! varname = 'agesno'
326	! agesno(:) = 0.0
327	! CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
328	! . jjm, rlonu, rlatv , interbar )
329
330	varname = 'zmea'
331	zmea(:) = 0.0
332	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
333	. jjm, rlonu, rlatv , interbar )
334
335	varname = 'zstd'
336	zstd(:) = 0.0
337	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
338	. jjm, rlonu, rlatv , interbar )
339	varname = 'zsig'
340	zsig(:) = 0.0
341	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
342	. jjm, rlonu, rlatv , interbar )
343	varname = 'zgam'
344	zgam(:) = 0.0
345	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
346	. jjm, rlonu, rlatv , interbar )
347	varname = 'zthe'
348	zthe(:) = 0.0
349	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
350	. jjm, rlonu, rlatv , interbar )
351	varname = 'zpic'
352	zpic(:) = 0.0
353	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
354	. jjm, rlonu, rlatv , interbar )
355	varname = 'zval'
356	zval(:) = 0.0
357	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
358	. jjm, rlonu, rlatv , interbar )
359	c
360	rugsrel(:) = 0.0
361	IF(ok_orodr) THEN
362	DO i = 1, iip1* jjp1
363	rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
364	ENDDO
365	ENDIF
366	C
367	C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
368	C utilise pour calculer les poids et pour assurer l'adequation entre les
369	C fractions d'ocean vu par l'atmosphere et l'ocean
370	C
371
372	write(,)'Essai de lecture masque ocean'
373	iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
374	if (iret .ne. 0) then
375	write(,)'ATTENTION!! pas de fichier o2a.nc trouve'
376	write(,)'Run force'
377	else
378	couple = .true.
379	iret = nf_close(nid_o2a)
380	call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
381	$ , nid_o2a)
382	if (iml_omask /= iim .or. jml_omask /= jjp1) then
383	write(,)'Dimensions non compatibles pour masque ocean'
384	write(,)'iim = ',iim,' iml_omask = ',iml_omask
385	write(,)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
386	stop
387	endif
388	ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
389	ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
390	ALLOCATE(dlon_omask(iml_omask), stat=iret)
391	ALLOCATE(dlat_omask(jml_omask), stat=iret)
392	ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
393	ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
394	CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
395	$ , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
396	CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
397	$ ttm_tmp, 1, 1, ocetmp)
398	CALL flinclo(fid)
399	dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
400	dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
401	ocemask = ocetmp
402	if (dlat_omask(1) < dlat_omask(jml_omask)) then
403	do j = 1, jml_omask
404	ocemask(:,j) = ocetmp(:,jml_omask-j+1)
405	enddo
406	endif
407	C
408	C passage masque ocean a la grille physique
409	C
410	ocemask_fi(1) = ocemask(1,1)
411	do j = 2, jjm
412	do i = 1, iim
413	ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
414	enddo
415	enddo
416	ocemask_fi(klon) = ocemask(1,jjp1)
417	zmasq = 1. - ocemask_fi
418	isst = 0
419	where (ocemask_fi(2:klon-1) >0.) isst = 1
420	endif
421
422
423	C
424	C lecture du fichier glace de terre pour fixer la fraction de terre
425	C et de glace de terre
426	C
427	CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
428	$ , fid)
429	ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
430	ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
431	ALLOCATE(dlon_lic(iml_lic), stat=iret)
432	ALLOCATE(dlat_lic(jml_lic), stat=iret)
433	ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
434	CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
435	$ , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
436	CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
437	$ , 1, 1, fraclic)
438	CALL flinclo(fid)
439	C
440	C interpolation sur la grille T du modele
441	C
442	WRITE(,) 'dimensions de landice iml_lic, jml_lic : ',
443	$ iml_lic, jml_lic
444	c
445	C sil les coordonnees sont en degres, on les transforme
446	C
447	IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) ) THEN
448	lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
449	ENDIF
450	IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
451	lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
452	ENDIF
453
454	dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
455	dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
456	C
457	CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
458	$ ,iim, jjp1,
459	$ rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
460	c$$$ flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
461	flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
462	C
463	C passage sur la grille physique
464	C
465	CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
466	$ pctsrf(1:klon, is_lic))
467	C adequation avec le maque terre/mer
468	c zmasq(157) = 0.
469	WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
470	pctsrf(1 : klon, is_lic) = 0.
471	END WHERE
472	WHERE (zmasq( 1 : klon) .LT. EPSFRA)
473	pctsrf(1 : klon, is_lic) = 0.
474	END WHERE
475	pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
476	DO ji = 1, klon
477	IF (zmasq(ji) .GT. EPSFRA) THEN
478	IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
479	pctsrf(ji, is_lic) = zmasq(ji)
480	pctsrf(ji, is_ter) = 0.
481	ELSE
482	pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
483	IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
484	pctsrf(ji,is_ter) = 0.
485	pctsrf(ji, is_lic) = zmasq(ji)
486	ENDIF
487	ENDIF
488	ENDIF
489	END DO
490	C
491	C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
492	C
493	pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
494
495
496	WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
497	pctsrf(1 : klon, is_oce) = 0.
498	END WHERE
499
500	if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
501
502	isst = 0
503	where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
504	C
505	C verif que somme des sous surface = 1
506	C
507	ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
508	$ .GT. EPSFRA)
509	IF (ji .NE. 0) THEN
510	WRITE(,) 'pb repartition sous maille pour ',ji,' points'
511	ENDIF
512
513	! where (pctsrf(1:klon, is_ter) >= .5)
514	! pctsrf(1:klon, is_ter) = 1.
515	! pctsrf(1:klon, is_oce) = 0.
516	! pctsrf(1:klon, is_sic) = 0.
517	! pctsrf(1:klon, is_lic) = 0.
518	! zmasq = 1.
519	! endwhere
520	! where (pctsrf(1:klon, is_lic) >= .5)
521	! pctsrf(1:klon, is_ter) = 0.
522	! pctsrf(1:klon, is_oce) = 0.
523	! pctsrf(1:klon, is_sic) = 0.
524	! pctsrf(1:klon, is_lic) = 1.
525	! zmasq = 1.
526	! endwhere
527	! where (pctsrf(1:klon, is_oce) >= .5)
528	! pctsrf(1:klon, is_ter) = 0.
529	! pctsrf(1:klon, is_oce) = 1.
530	! pctsrf(1:klon, is_sic) = 0.
531	! pctsrf(1:klon, is_lic) = 0.
532	! zmasq = 0.
533	! endwhere
534	! where (pctsrf(1:klon, is_sic) >= .5)
535	! pctsrf(1:klon, is_ter) = 0.
536	! pctsrf(1:klon, is_oce) = 0.
537	! pctsrf(1:klon, is_sic) = 1.
538	! pctsrf(1:klon, is_lic) = 0.
539	! zmasq = 0.
540	! endwhere
541	! call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
542	C
543	C verif que somme des sous surface = 1
544	C
545	! ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
546	! $ .GT. EPSFRA)
547	! IF (ji .NE. 0) THEN
548	! WRITE(,) 'pb repartition sous maille pour ',ji,' points'
549	! ENDIF
550
551
552
553
554	C Calcul intermediaire
555	c
556	CALL massdair( p3d, masse )
557	c
558
559	print *,' ALPHAX ',alphax
560
561	DO l = 1, llm
562	DO i = 1, iim
563	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
564	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
565	ENDDO
566	xpn = SUM(xppn)/apoln
567	xps = SUM(xpps)/apols
568	DO i = 1, iip1
569	masse( i , 1 , l ) = xpn
570	masse( i , jjp1 , l ) = xps
571	ENDDO
572	ENDDO
573	q3d(iip1,:,:,:) = q3d(1,:,:,:)
574	phis(iip1,:) = phis(1,:)
575
576	C Ecriture
577
578
579	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
580	* tetagdiv, tetagrot , tetatemp )
581	print*,'sortie inidissip'
582	itau = 0
583	itau_dyn = 0
584	itau_phy = 0
585	iday = dayref +itau/day_step
586	time = FLOAT(itau-(iday-dayref)*day_step)/day_step
587	c
588	IF(time.GT.1) THEN
589	time = time - 1
590	iday = iday + 1
591	ENDIF
592	day_ref = dayref
593	annee_ref = anneeref
594
595	CALL geopot ( ip1jmp1, tpot , pk , pks, phis , phi )
596	print*,'sortie geopot'
597
598	CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
599	* phi,w, pbaru,pbarv,time+iday-dayref )
600	print*,'sortie caldyn0'
601	CALL dynredem0("start.nc",dayref,phis,nqmx)
602	print*,'sortie dynredem0'
603	CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
604	. psol)
605	print*,'sortie dynredem1'
606	C
607	C Ecriture etat initial physique
608	C
609	phystep = dtvr * FLOAT(iphysiq)
610	radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) )
611	co2_ppm = 330.0
612	solaire = 1370.0
613
614	c
615	c Initialisation
616	c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
617	c
618	tsolsrf(:,is_ter) = tsol
619	tsolsrf(:,is_lic) = tsol
620	tsolsrf(:,is_oce) = tsol
621	tsolsrf(:,is_sic) = tsol
622	snsrf(:,is_ter) = sn
623	snsrf(:,is_lic) = sn
624	snsrf(:,is_oce) = sn
625	snsrf(:,is_sic) = sn
626	albe(:,is_ter) = 0.08
627	albe(:,is_lic) = 0.6
628	albe(:,is_oce) = 0.5
629	albe(:,is_sic) = 0.6
630	evap(:,:) = 0.
631	qsolsrf(:,is_ter) = 150
632	qsolsrf(:,is_lic) = 150
633	qsolsrf(:,is_oce) = 150.
634	qsolsrf(:,is_sic) = 150.
635	do i = 1, nbsrf
636	do j = 1, nsoilmx
637	tsoil(:,j,i) = tsol
638	enddo
639	enddo
640	rain_fall = 0.; snow_fall = 0.
641	solsw = 165.
642	sollw = -53.
643	t_ancien = 273.15
644	q_ancien = 0.
645	agesno = 0.
646	deltat = 0.
647	frugs(1:klon,is_oce) = rugmer(1:klon)
648	frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
649	frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
650	frugs(1:klon,is_sic) = 0.001
651	fder = 0.0
652
653	call phyredem("startphy.nc",phystep,radpas, co2_ppm, solaire,
654	$ latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, snsrf,
655	$ albe, evap, rain_fall, snow_fall, solsw, sollw, fder,
656	$ radsol, frugs, agesno,
657	$ zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
658	$ t_ancien, q_ancien)
659
660	C Sortie Visu pour les champs dynamiques
661	print*,'sortie visu'
662	time_step = 1.
663	t_ops = 2.
664	t_wrt = 2.
665	itau = 2.
666	visu_file='Etat0_visu.nc'
667	CALL initdynav(visu_file,dayref,anneeref,time_step,
668	. t_ops, t_wrt, nqmx, visuid)
669	CALL writedynav(visuid, nqmx, itau,vvent ,
670	. uvent,tpot,pk,phi,q3d,masse,psol,phis)
671	print*,'entree histclo'
672	CALL histclo
673	RETURN
674	!
675	END SUBROUTINE etat0_netcdf

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

Download in other formats:

Original Format