Context Navigation

etat0_netcdf.F @ 1132

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

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

Download in other formats:

Original Format