Context Navigation

etat0_netcdf.F @ 1116

Last change on this file since 1116 was 1114, checked in by jghattas, 17 years ago

Creation du module infotrac:

contient les variables de advtrac.h
contient la subroutine iniadvtrac renommer en infotrac_init
le nombre des traceurs est lu dans tracer.def en dynamique (ou par default ou recu par INCA)
ce module est utilise dans la dynamique et la physique
contient aussi la variable nbtr qui avant etait stockee dans dimphy

Le fichier advtrac.h n'existe plus.
La compilation ne prend plus en compte le nombre de traceur.

/JG

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 23.3 KB

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

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

Download in other formats:

Original Format