Context Navigation

etat0_netcdf.F @ 1096

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

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

Download in other formats:

Original Format