Context Navigation

etat0_netcdf.F @ 1198

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

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

Download in other formats:

Original Format