Context Navigation

etat0_netcdf.F @ 1033

Last change on this file since 1033 was 995, checked in by lsce, 16 years ago

Modifications liées au calcul des nouveau sous-fractions
Nettoyage de ocean slab : il reste uniquement la version avec glace de mer forcé
Nouveaux variables pour distiguer la version et type d'ocean : type_ocean=force/slab/couple, version_ocean=opa8/nemo pour couplé ou version_ocean=sicOBS pour slab

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats:

Original Format