Context Navigation

etat0_netcdf.F @ 5133

Last change on this file since 5133 was 546, checked in by lmdzadmin, 20 years ago
Suite des adaptations au portage sur SGI et VPP pour Prism (Caubel & Demory) LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.1 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), qsurf(klon,nbsrf),snsrf(klon,nbsrf)
38	REAL :: albe(klon,nbsrf), evap(klon,nbsrf), alblw(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), run_off_lic_0(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	real :: clwcon(klon,klev),rnebcon(klon,klev),ratqs(klon,klev)
51	! declarations pour lecture glace de mer
52	REAL :: rugv(klon)
53	INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
54	INTEGER :: itaul(1), fid
55	REAL :: lev(1), date
56	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
57	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic
58	REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
59	REAL :: flic_tmp(iip1, jjp1)
60	REAL :: champint(iim, jjp1)
61	!
62
63	CHARACTER*80 :: varname
64	!
65	INTEGER :: i,j, ig, l, ji,ii1,ii2, n
66	REAL :: xpi
67	!
68	REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
69	REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
70	REAL :: workvar(iip1,jjp1,llm)
71	!
72	REAL :: prefkap, unskap
73	!
74	REAL :: q_sat
75	EXTERNAL q_sat
76	real :: time_step,t_ops,t_wrt
77
78	#include "comdissnew.h"
79	#include "control.h"
80	#include "serre.h"
81	#include "clesphys.h"
82
83	INTEGER :: longcles
84	PARAMETER ( longcles = 20 )
85	REAL :: clesphy0 ( longcles )
86	REAL :: p(iip1,jjp1,llm)
87	INTEGER :: itau, iday
88	REAL :: masse(iip1,jjp1,llm)
89	REAL :: xpn,xps,xppn(iim),xpps(iim)
90	real :: time
91	REAL :: phi(ip1jmp1,llm)
92	REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
93	REAL :: w(ip1jmp1,llm)
94	cIM REAL ::phystep,co2_ppm,solaire
95	REAL ::phystep
96	INTEGER :: radpas
97	real zrel(iip1*jjp1),chmin,chmax
98
99	CHARACTER*80 :: visu_file
100	INTEGER :: visuid
101
102	! pour la lecture du fichier masque ocean
103	integer :: nid_o2a
104	logical :: couple = .false.
105	INTEGER :: iml_omask, jml_omask
106	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
107	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask
108	REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
109	real, dimension(klon) :: ocemask_fi
110	integer :: isst(klon-2)
111	real zx_tmp_2d(iim,jjp1)
112	!
113	! Constantes
114	!
115	pi = 4. * ATAN(1.)
116	rad = 6371229.
117	omeg = 4.* ASIN(1.)/(24.*3600.)
118	g = 9.8
119	daysec = 86400.
120	kappa = 0.2857143
121	cpp = 1004.70885
122	!
123	preff = 101325.
124	unskap = 1./kappa
125	!
126	jmp1 = jjm + 1
127	!
128	! Construct a grid
129	!
130
131	! CALL defrun_new(99,.TRUE.,clesphy0)
132	CALL conf_gcm( 99, .TRUE. , clesphy0 )
133
134	dtvr = daysec/FLOAT(day_step)
135	print*,'dtvr',dtvr
136
137	CALL inicons0()
138	CALL inigeom()
139	!
140	CALL inifilr()
141	!
142	latfi(1) = ASIN(1.0)
143	DO j = 2, jjm
144	DO i = 1, iim
145	latfi((j-2)*iim+1+i)= rlatu(j)
146	ENDDO
147	ENDDO
148	latfi(klon) = - ASIN(1.0)
149	!
150	lonfi(1) = 0.0
151	DO j = 2, jjm
152	DO i = 1, iim
153	lonfi((j-2)*iim+1+i) = rlonv(i)
154	ENDDO
155	ENDDO
156	lonfi(klon) = 0.0
157	!
158	xpi = 2.0 * ASIN(1.0)
159	DO ig = 1, klon
160	latfi(ig) = latfi(ig) * 180.0 / xpi
161	lonfi(ig) = lonfi(ig) * 180.0 / xpi
162	ENDDO
163	!
164
165
166	C
167	C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
168	C utilise pour calculer les poids et pour assurer l'adequation entre les
169	C fractions d'ocean vu par l'atmosphere et l'ocean. Sinon, on cree le masque
170	C a partir du fichier relief
171	C
172
173	write(,)'Essai de lecture masque ocean'
174	iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
175	if (iret .ne. 0) then
176	write(,)'ATTENTION!! pas de fichier o2a.nc trouve'
177	write(,)'Run force'
178	varname = 'masque'
179	masque(:,:) = 0.0
180	CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
181	, jjm ,rlonu,rlatv , interbar )
182	WRITE(,) 'MASQUE construit : Masque'
183	WRITE(*,'(97I1)') nINT(masque(:,:))
184	call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
185	WHERE (zmasq(1 : klon) .LT. EPSFRA)
186	zmasq(1 : klon) = 0.
187	END WHERE
188	WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
189	zmasq(1 : klon) = 1.
190	END WHERE
191	else
192	couple = .true.
193	iret = nf_close(nid_o2a)
194	call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
195	$ , nid_o2a)
196	if (iml_omask /= iim .or. jml_omask /= jjp1) then
197	write(,)'Dimensions non compatibles pour masque ocean'
198	write(,)'iim = ',iim,' iml_omask = ',iml_omask
199	write(,)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
200	stop
201	endif
202	ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
203	ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
204	ALLOCATE(dlon_omask(iml_omask), stat=iret)
205	ALLOCATE(dlat_omask(jml_omask), stat=iret)
206	ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
207	ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
208	CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
209	$ , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
210	CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
211	$ ttm_tmp, 1, 1, ocetmp)
212	CALL flinclo(fid)
213	dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
214	dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
215	ocemask = ocetmp
216	if (dlat_omask(1) < dlat_omask(jml_omask)) then
217	do j = 1, jml_omask
218	ocemask(:,j) = ocetmp(:,jml_omask-j+1)
219	enddo
220	endif
221	C
222	C passage masque ocean a la grille physique
223	C
224	write(,)'ocemask '
225	write(*,'(96i1)')int(ocemask)
226	ocemask_fi(1) = ocemask(1,1)
227	do j = 2, jjm
228	do i = 1, iim
229	ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
230	enddo
231	enddo
232	ocemask_fi(klon) = ocemask(1,jjp1)
233	zmasq = 1. - ocemask_fi
234	endif
235
236	call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
237
238	varname = 'relief'
239	! This line needs to be replaced by a call to restget to get the values in the restart file
240	orog(:,:) = 0.0
241	CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
242	, jjm ,rlonu,rlatv , interbar, masque )
243	!
244	WRITE(,) 'OUT OF GET VARIABLE : Relief'
245	! WRITE(*,'(49I1)') INT(orog(:,:))
246	!
247	varname = 'rugosite'
248	! This line needs to be replaced by a call to restget to get the values in the restart file
249	rugo(:,:) = 0.0
250	CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
251	, jjm, rlonu,rlatv , interbar )
252	!
253	WRITE(,) 'OUT OF GET VARIABLE : Rugosite'
254	! WRITE(,'(49I1)') INT(rugo(:,:)10)
255	!
256	C
257	C on initialise les sous surfaces
258	C
259	pctsrf=0.
260	c
261	varname = 'psol'
262	psol(:,:) = 0.0
263	CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
264	, jjm ,rlonu,rlatv , interbar )
265	!
266	! Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
267	! anyway.
268	!
269	! WRITE(,) 'PSOL :', psol(10,20)
270	! WRITE(,) ap(:), bp(:)
271	CALL pression(ip1jmp1, ap, bp, psol, p3d)
272	! WRITE(,) 'P3D :', p3d(10,20,:)
273	CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
274	! WRITE(,) 'PK:', pk(10,20,:)
275	!
276	!
277	!
278	prefkap = preff ** kappa
279	! WRITE(,) 'unskap, cpp, preff :', unskap, cpp, preff
280	DO l = 1, llm
281	DO j=1,jjp1
282	DO i =1, iip1
283	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
284	ENDDO
285	ENDDO
286	ENDDO
287	!
288	! WRITE(,) 'PLS :', pls(10,20,:)
289	!
290	varname = 'surfgeo'
291	phis(:,:) = 0.0
292	CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
293	, jjm ,rlonu,rlatv, interbar )
294	!
295	varname = 'u'
296	uvent(:,:,:) = 0.0
297	CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
298	. workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
299	!
300	varname = 'v'
301	vvent(:,:,:) = 0.0
302	CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
303	. workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
304	!
305	varname = 't'
306	t3d(:,:,:) = 0.0
307	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
308	. workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
309	!
310	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
311	. maxval(t3d(:,:,:))
312	varname = 'tpot'
313	tpot(:,:,:) = 0.0
314	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
315	. pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
316	!
317	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
318	. maxval(t3d(:,:,:))
319	WRITE(,) 'PLS min,max:',minval(pls(:,:,:)),
320	. maxval(pls(:,:,:))
321	DO l = 1, llm
322	DO j=1,jjp1
323	DO i =1, iip1-1
324	qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
325	ENDDO
326	qsat(iip1,j,l) = qsat(1,j,l)
327	ENDDO
328	ENDDO
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 Ecriture
580
581
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.77
634	albe(:,is_oce) = 0.5
635	albe(:,is_sic) = 0.6
636	alblw(:,is_ter) = albe(:,is_ter)
637	alblw(:,is_lic) = albe(:,is_lic)
638	alblw(:,is_oce) = albe(:,is_oce)
639	alblw(:,is_sic) = albe(:,is_sic)
640
641	evap(:,:) = 0.
642	C vapeur d'eau a la surface =- vapeur d'eau premier niveau
643	do n = 1, nbsrf
644	DO j=2,jjm
645	DO i =1, iim
646	qsurf((j-2)*iim+1+i,n)=q3d(i,j,1,1)
647	END DO
648	END DO
649	qsurf(1,n)=q3d(1,1,1,1) ! pole nord
650	qsurf(klon,n)=q3d(1,jjp1,1,1) ! pole sud
651	END DO
652	qsol(:) = 150.
653	do i = 1, nbsrf
654	do j = 1, nsoilmx
655	tsoil(:,j,i) = tsol
656	enddo
657	enddo
658	rain_fall = 0.; snow_fall = 0.
659	solsw = 165.
660	sollw = -53.
661	t_ancien = 273.15
662	q_ancien = 0.
663	agesno = 0.
664	deltat = 0.
665	frugs(1:klon,is_oce) = rugmer(1:klon)
666	frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
667	frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
668	frugs(1:klon,is_sic) = 0.002
669	fder = 0.0
670	clwcon = 0.0
671	rnebcon = 0.0
672	ratqs = 0.0
673	run_off_lic_0 = 0.0
674
675	cIM call phyredem("startphy.nc",phystep,radpas, co2_ppm, solaire,
676	call phyredem("startphy.nc",phystep,radpas,
677	$ latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsurf, qsol,
678	$ snsrf, albe, alblw,
679	$ evap, rain_fall, snow_fall, solsw, sollw, fder,
680	$ radsol, frugs, agesno,
681	$ zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
682	$ t_ancien, q_ancien, rnebcon, ratqs, clwcon,
683	$ run_off_lic_0)
684
685	C Sortie Visu pour les champs dynamiques
686	print*,'sortie visu'
687	time_step = 1.
688	t_ops = 2.
689	t_wrt = 2.
690	itau = 2.
691	visu_file='Etat0_visu.nc'
692	CALL initdynav(visu_file,dayref,anneeref,time_step,
693	. t_ops, t_wrt, nqmx, visuid)
694	CALL writedynav(visuid, nqmx, itau,vvent ,
695	. uvent,tpot,pk,phi,q3d,masse,psol,phis)
696	print*,'entree histclo'
697	CALL histclo
698	RETURN
699	!
700	END SUBROUTINE etat0_netcdf

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

Download in other formats:

Original Format