Context Navigation

etat0_netcdf.F @ 1261

Last change on this file since 1261 was 1245, checked in by Laurent Fairhead, 15 years ago
Un oubli dans le conf_phys suite a la correction de bug sur ratqs
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 24.2 KB

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

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

Download in other formats:

Original Format