Context Navigation

etat0_netcdf.F @ 1295

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

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

Download in other formats:

Original Format