Context Navigation

etat0_netcdf.F @ 988

Last change on this file since 988 was 988, checked in by Laurent Fairhead, 16 years ago
Corrections suite a la re-ecriture de phyredem.F Quelques initialisations de plus LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 23.9 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 ocean_slab_mod
14	USE surface_data, ONLY : ocean, ok_veget
15	USE phys_state_var_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	. 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,nqmx), 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	cIM "slab" ocean
49	REAL :: tslab(klon), seaice(klon)
50	REAL :: frugs(klon,nbsrf), agesno(klon,nbsrf)
51	REAL :: rugmer(klon)
52	REAL :: qd(iip1, jjp1, llm)
53	REAL :: run_off_lic_0(klon)
54	! declarations pour lecture glace de mer
55	REAL :: rugv(klon)
56	INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
57	INTEGER :: itaul(1), fid
58	REAL :: lev(1), date
59	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
60	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic
61	REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
62	REAL :: flic_tmp(iip1, jjp1)
63	REAL :: champint(iim, jjp1)
64	!
65
66	CHARACTER*80 :: varname
67	!
68	INTEGER :: i,j, ig, l, ji,ii1,ii2
69	INTEGER :: nq
70	REAL :: xpi
71	!
72	REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
73	REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
74	REAL :: workvar(iip1,jjp1,llm)
75	!
76	REAL :: prefkap, unskap
77	!
78	real :: time_step,t_ops,t_wrt
79
80	#include "comdissnew.h"
81	#include "control.h"
82	#include "serre.h"
83	#include "clesphys.h"
84
85	INTEGER :: longcles
86	PARAMETER ( longcles = 20 )
87	REAL :: clesphy0 ( longcles )
88	REAL :: p(iip1,jjp1,llm)
89	INTEGER :: itau, iday
90	REAL :: masse(iip1,jjp1,llm)
91	REAL :: xpn,xps,xppn(iim),xpps(iim)
92	real :: time
93	REAL :: phi(ip1jmp1,llm)
94	REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
95	REAL :: w(ip1jmp1,llm)
96	REAL ::phystep
97	REAL :: rugsrel(iip1*jjp1)
98	REAL :: fder(klon)
99	real zrel(iip1*jjp1),chmin,chmax
100
101	CHARACTER*80 :: visu_file
102	INTEGER :: visuid
103
104	! pour la lecture du fichier masque ocean
105	integer :: nid_o2a
106	logical :: couple = .false.
107	INTEGER :: iml_omask, jml_omask
108	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
109	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask
110	REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
111	real, dimension(klon) :: ocemask_fi
112	integer :: isst(klon-2)
113	real zx_tmp_2d(iim,jjp1)
114
115	REAL :: dummy(klon, nbsrf)
116
117	logical :: ok_newmicro
118	integer :: iflag_radia
119	logical :: ok_journe, ok_mensuel, ok_instan, ok_hf
120	LOGICAL :: ok_ade, ok_aie, aerosol_couple
121	REAL :: bl95_b0, bl95_b1
122	real :: fact_cldcon, facttemps,ratqsbas,ratqshaut
123	integer :: iflag_cldcon
124	integer :: iflag_ratqs
125	integer :: iflag_coupl
126	integer :: iflag_clos
127	integer :: iflag_wake
128	integer :: iflag_thermals,nsplit_thermals
129	real :: tau_thermals
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(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &
155	& ok_hf, &
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_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	cIM "slab" ocean
421	varname = 'tslab'
422	tslab(:) = 0.0
423	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,tslab,0.0,
424	. jjm, rlonu, rlatv , interbar )
425	c
426	varname = 'seaice'
427	seaice(:) = 0.0
428	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,seaice,0.0,
429	. jjm, rlonu, rlatv , interbar )
430	!
431	varname = 'rugmer'
432	rugmer(:) = 0.0
433	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
434	. jjm, rlonu, rlatv , interbar )
435	!
436	! varname = 'agesno'
437	! agesno(:) = 0.0
438	! CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
439	! . jjm, rlonu, rlatv , interbar )
440
441	varname = 'zmea'
442	zmea(:) = 0.0
443	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
444	. jjm, rlonu, rlatv , interbar )
445
446	varname = 'zstd'
447	zstd(:) = 0.0
448	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
449	. jjm, rlonu, rlatv , interbar )
450	varname = 'zsig'
451	zsig(:) = 0.0
452	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
453	. jjm, rlonu, rlatv , interbar )
454	varname = 'zgam'
455	zgam(:) = 0.0
456	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
457	. jjm, rlonu, rlatv , interbar )
458	varname = 'zthe'
459	zthe(:) = 0.0
460	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
461	. jjm, rlonu, rlatv , interbar )
462	varname = 'zpic'
463	zpic(:) = 0.0
464	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
465	. jjm, rlonu, rlatv , interbar )
466	varname = 'zval'
467	zval(:) = 0.0
468	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
469	. jjm, rlonu, rlatv , interbar )
470	c
471	rugsrel(:) = 0.0
472	IF(ok_orodr) THEN
473	DO i = 1, iip1* jjp1
474	rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
475	ENDDO
476	ENDIF
477
478
479	C
480	C lecture du fichier glace de terre pour fixer la fraction de terre
481	C et de glace de terre
482	C
483	CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
484	$ , fid)
485	ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
486	ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
487	ALLOCATE(dlon_lic(iml_lic), stat=iret)
488	ALLOCATE(dlat_lic(jml_lic), stat=iret)
489	ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
490	CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
491	$ , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
492	CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
493	$ , 1, 1, fraclic)
494	CALL flinclo(fid)
495	C
496	C interpolation sur la grille T du modele
497	C
498	WRITE(,) 'dimensions de landice iml_lic, jml_lic : ',
499	$ iml_lic, jml_lic
500	c
501	C sil les coordonnees sont en degres, on les transforme
502	C
503	IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) ) THEN
504	lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
505	ENDIF
506	IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
507	lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
508	ENDIF
509
510	dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
511	dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
512	C
513	CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
514	$ ,iim, jjp1,
515	$ rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
516	cx$$$ flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
517	flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
518	C
519	C passage sur la grille physique
520	C
521	CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
522	$ pctsrf(1:klon, is_lic))
523	C adequation avec le maque terre/mer
524	c zmasq(157) = 0.
525	WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
526	pctsrf(1 : klon, is_lic) = 0.
527	END WHERE
528	WHERE (zmasq( 1 : klon) .LT. EPSFRA)
529	pctsrf(1 : klon, is_lic) = 0.
530	END WHERE
531	pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
532	DO ji = 1, klon
533	IF (zmasq(ji) .GT. EPSFRA) THEN
534	IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
535	pctsrf(ji, is_lic) = zmasq(ji)
536	pctsrf(ji, is_ter) = 0.
537	ELSE
538	pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
539	IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
540	pctsrf(ji,is_ter) = 0.
541	pctsrf(ji, is_lic) = zmasq(ji)
542	ENDIF
543	ENDIF
544	ENDIF
545	END DO
546	C
547	C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
548	C
549	pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
550
551
552	WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
553	pctsrf(1 : klon, is_oce) = 0.
554	END WHERE
555
556	if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
557
558	isst = 0
559	where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
560	C
561	C verif que somme des sous surface = 1
562	C
563	ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf),dim=2))-1.0)
564	$ .GT. EPSFRA)
565	IF (ji .NE. 0) THEN
566	WRITE(,) 'pb repartition sous maille pour ',ji,' points'
567	ENDIF
568
569	! where (pctsrf(1:klon, is_ter) >= .5)
570	! pctsrf(1:klon, is_ter) = 1.
571	! pctsrf(1:klon, is_oce) = 0.
572	! pctsrf(1:klon, is_sic) = 0.
573	! pctsrf(1:klon, is_lic) = 0.
574	! zmasq = 1.
575	! endwhere
576	! where (pctsrf(1:klon, is_lic) >= .5)
577	! pctsrf(1:klon, is_ter) = 0.
578	! pctsrf(1:klon, is_oce) = 0.
579	! pctsrf(1:klon, is_sic) = 0.
580	! pctsrf(1:klon, is_lic) = 1.
581	! zmasq = 1.
582	! endwhere
583	! where (pctsrf(1:klon, is_oce) >= .5)
584	! pctsrf(1:klon, is_ter) = 0.
585	! pctsrf(1:klon, is_oce) = 1.
586	! pctsrf(1:klon, is_sic) = 0.
587	! pctsrf(1:klon, is_lic) = 0.
588	! zmasq = 0.
589	! endwhere
590	! where (pctsrf(1:klon, is_sic) >= .5)
591	! pctsrf(1:klon, is_ter) = 0.
592	! pctsrf(1:klon, is_oce) = 0.
593	! pctsrf(1:klon, is_sic) = 1.
594	! pctsrf(1:klon, is_lic) = 0.
595	! zmasq = 0.
596	! endwhere
597	! call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
598	C
599	C verif que somme des sous surface = 1
600	C
601	! ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
602	! $ .GT. EPSFRA)
603	! IF (ji .NE. 0) THEN
604	! WRITE(,) 'pb repartition sous maille pour ',ji,' points'
605	! ENDIF
606
607	CALL gr_fi_ecrit(1,klon,iim,jjp1,zmasq,zx_tmp_2d)
608	write(,)'zmasq = '
609	write(*,'(96i1)')nint(zx_tmp_2d)
610	call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
611	WRITE(,) 'MASQUE construit : Masque'
612	WRITE(*,'(97I1)') nINT(masque(:,:))
613
614
615
616	C Calcul intermediaire
617	c
618	CALL massdair( p3d, masse )
619	c
620
621	print *,' ALPHAX ',alphax
622
623	DO l = 1, llm
624	DO i = 1, iim
625	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
626	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
627	ENDDO
628	xpn = SUM(xppn)/apoln
629	xps = SUM(xpps)/apols
630	DO i = 1, iip1
631	masse( i , 1 , l ) = xpn
632	masse( i , jjp1 , l ) = xps
633	ENDDO
634	ENDDO
635	q3d(iip1,:,:,:) = q3d(1,:,:,:)
636	phis(iip1,:) = phis(1,:)
637
638	C init pour traceurs
639	call iniadvtrac(nq)
640	C Ecriture
641	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
642	* tetagdiv, tetagrot , tetatemp )
643	print*,'sortie inidissip'
644	itau = 0
645	itau_dyn = 0
646	itau_phy = 0
647	iday = dayref +itau/day_step
648	time = FLOAT(itau-(iday-dayref)*day_step)/day_step
649	c
650	IF(time.GT.1) THEN
651	time = time - 1
652	iday = iday + 1
653	ENDIF
654	day_ref = dayref
655	annee_ref = anneeref
656
657	CALL geopot ( ip1jmp1, tpot , pk , pks, phis , phi )
658	print*,'sortie geopot'
659
660	CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
661	* phi,w, pbaru,pbarv,time+iday-dayref )
662	print*,'sortie caldyn0'
663	CALL dynredem0("start.nc",dayref,phis,nqmx)
664	print*,'sortie dynredem0'
665	CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
666	. psol)
667	print*,'sortie dynredem1'
668	C
669	C Ecriture etat initial physique
670	C
671	write(,)'phystep ',dtvr,iphysiq,nbapp_rad
672	phystep = dtvr * FLOAT(iphysiq)
673	radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) )
674	write(,)'phystep =', phystep, radpas
675	cIM : lecture de co2_ppm & solaire ds physiq.def
676	c co2_ppm = 348.0
677	c solaire = 1365.0
678
679	c
680	c Initialisation
681	c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
682	c
683	ftsol(:,is_ter) = tsol
684	ftsol(:,is_lic) = tsol
685	ftsol(:,is_oce) = tsol
686	ftsol(:,is_sic) = tsol
687	snsrf(:,is_ter) = sn
688	snsrf(:,is_lic) = sn
689	snsrf(:,is_oce) = sn
690	snsrf(:,is_sic) = sn
691	falb1(:,is_ter) = 0.08
692	falb1(:,is_lic) = 0.6
693	falb1(:,is_oce) = 0.5
694	falb1(:,is_sic) = 0.6
695	falb2 = falb1
696	evap(:,:) = 0.
697	qsolsrf(:,is_ter) = 150
698	qsolsrf(:,is_lic) = 150
699	qsolsrf(:,is_oce) = 150.
700	qsolsrf(:,is_sic) = 150.
701	do i = 1, nbsrf
702	do j = 1, nsoilmx
703	tsoil(:,j,i) = tsol
704	enddo
705	enddo
706	rain_fall = 0.; snow_fall = 0.
707	solsw = 165.
708	sollw = -53.
709	t_ancien = 273.15
710	q_ancien = 0.
711	agesno = 0.
712	cIM "slab" ocean
713	tslab(1:klon) = ftsol(1:klon,is_oce)
714	seaice = 0. !kg m-2
715	c
716	frugs(1:klon,is_oce) = rugmer(1:klon)
717	frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
718	frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
719	frugs(1:klon,is_sic) = 0.001
720	fder = 0.0
721	clwcon = 0.0
722	rnebcon = 0.0
723	ratqs = 0.0
724	run_off_lic_0 = 0.0
725	rugoro = 0.0
726
727	c
728	c Avant l'appel a phyredem, on initialize les modules de surface
729	c avec les valeurs qui vont etre ecrit dans startphy.nc
730	c
731	dummy(:,:) = 1.0
732	ocean = "slab"
733	pbl_tke(:,:,:) = 1.e-8
734	zmax0(:) = 40.
735	f0(:) = 1.e-5
736	ema_work1(:,:) = 0.
737	ema_work2(:,:) = 0.
738	wake_deltat(:,:) = 0.
739	wake_deltaq(:,:) = 0.
740	wake_s(:) = 0.
741	wake_cstar(:) = 0.
742	wake_fip(:) = 0.
743
744	call fonte_neige_init(run_off_lic_0)
745	call pbl_surface_init(qsol, fder, snsrf, qsolsrf,
746	$ evap, frugs, agesno, tsoil)
747	call ocean_slab_init(dummy(1,1), tslab, seaice, dummy(:,:))
748
749	call phyredem("startphy.nc")
750
751
752
753	C Sortie Visu pour les champs dynamiques
754	if (1.eq.0 ) then
755	print*,'sortie visu'
756	time_step = 1.
757	t_ops = 2.
758	t_wrt = 2.
759	itau = 2.
760	visu_file='Etat0_visu.nc'
761	CALL initdynav(visu_file,dayref,anneeref,time_step,
762	. t_ops, t_wrt, nqmx, visuid)
763	CALL writedynav(visuid, nqmx, itau,vvent ,
764	. uvent,tpot,pk,phi,q3d,masse,psol,phis)
765	else
766	print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
767	endif
768	print*,'entree histclo'
769	CALL histclo
770	RETURN
771	!
772	END SUBROUTINE etat0_netcdf
773

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

Download in other formats:

Original Format