Context Navigation

etat0_netcdf.F @ 1119

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

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

Download in other formats:

Original Format