Context Navigation

etat0_netcdf.F90 @ 4057

Last change on this file since 4057 was 1492, checked in by Laurent Fairhead, 14 years ago

Merge of development branch LMDZ5V2.0-dev r1455:r1491 into the trunk.
Validation made locally: restart files are strictly equal between the HEAD of the trunk
and r1491 of LMDZ5V2.0-dev

Synchro de la branche de développement LMDZ5V2.0-dev r1455:r1491 et de la trunk
Validation faite en local: les fichiers restart sont équivalents entre la HEAD de la trunk
et la révision r1491 de LMDZ5V2.0-dev

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

File size: 19.3 KB

Line
1	!
2	! $Id: etat0_netcdf.F90 1492 2011-03-08 08:10:25Z dcugnet $
3	!
4	!-------------------------------------------------------------------------------
5	!
6	SUBROUTINE etat0_netcdf(ib, masque, letat0)
7	!
8	!-------------------------------------------------------------------------------
9	! Purpose: Creates initial states
10	!-------------------------------------------------------------------------------
11	! Note: This routine is designed to work for Earth
12	!-------------------------------------------------------------------------------
13	USE control_mod
14	#ifdef CPP_EARTH
15	USE startvar
16	USE ioipsl
17	USE dimphy
18	USE infotrac
19	USE fonte_neige_mod
20	USE pbl_surface_mod
21	USE phys_state_var_mod
22	USE filtreg_mod
23	USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
24	USE conf_phys_m, ONLY: conf_phys
25	! For parameterization of ozone chemistry:
26	use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
27	use press_coefoz_m, only: press_coefoz
28	use regr_pr_o3_m, only: regr_pr_o3
29	USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
30	#endif
31	IMPLICIT NONE
32	!-------------------------------------------------------------------------------
33	! Arguments:
34	#include "dimensions.h"
35	#include "paramet.h"
36	#include "iniprint.h"
37	LOGICAL, INTENT(IN) :: ib ! barycentric interpolat.
38	REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
39	LOGICAL, INTENT(IN) :: letat0 ! F: masque only required
40	#ifndef CPP_EARTH
41	WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
42	#else
43	!-------------------------------------------------------------------------------
44	! Local variables:
45	#include "comgeom2.h"
46	#include "comvert.h"
47	#include "comconst.h"
48	#include "indicesol.h"
49	#include "dimsoil.h"
50	#include "temps.h"
51	REAL, DIMENSION(klon) :: tsol, qsol
52	REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0
53	REAL, DIMENSION(iip1,jjp1) :: orog, rugo, psol, phis
54	REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
55	REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd
56	REAL, DIMENSION(iip1,jjm ,llm) :: vvent
57	REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: q3d
58	REAL, DIMENSION(klon,nbsrf) :: qsolsrf, snsrf, evap
59	REAL, DIMENSION(klon,nbsrf) :: frugs, agesno
60	REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
61
62	!--- Local variables for sea-ice reading:
63	INTEGER :: iml_lic, jml_lic, llm_tmp
64	INTEGER :: ttm_tmp, iret, fid
65	INTEGER, DIMENSION(1) :: itaul
66	REAL, DIMENSION(1) :: lev
67	REAL :: date
68	REAL, DIMENSION(:,:), ALLOCATABLE :: lon_lic, lat_lic, fraclic
69	REAL, DIMENSION(:), ALLOCATABLE :: dlon_lic, dlat_lic
70	REAL, DIMENSION(iip1,jjp1) :: flic_tmp
71
72	!--- Misc
73	CHARACTER(LEN=80) :: x, fmt
74	INTEGER :: i, j, l, ji
75	REAL, DIMENSION(iip1,jjp1,llm) :: alpha, beta, pk, pls, y
76	REAL, DIMENSION(ip1jmp1) :: pks
77
78	#include "comdissnew.h"
79	#include "serre.h"
80	#include "clesphys.h"
81
82	REAL, DIMENSION(iip1,jjp1,llm) :: masse
83	INTEGER :: itau, iday
84	REAL :: xpn, xps, time, phystep
85	REAL, DIMENSION(iim) :: xppn, xpps
86	REAL, DIMENSION(ip1jmp1,llm) :: pbaru, phi, w
87	REAL, DIMENSION(ip1jm ,llm) :: pbarv
88	REAL, DIMENSION(klon) :: fder
89
90	!--- Local variables for ocean mask reading:
91	INTEGER :: nid_o2a, iml_omask, jml_omask
92	LOGICAL :: couple=.FALSE.
93	REAL, DIMENSION(:,:), ALLOCATABLE :: lon_omask, lat_omask, ocemask, ocetmp
94	REAL, DIMENSION(:), ALLOCATABLE :: dlon_omask,dlat_omask
95	REAL, DIMENSION(klon) :: ocemask_fi
96	INTEGER, DIMENSION(klon-2) :: isst
97	REAL, DIMENSION(iim,jjp1) :: zx_tmp_2d
98	REAL :: dummy
99	LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf
100	LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats
101	INTEGER :: iflag_radia, flag_aerosol
102	REAL :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut
103	REAL :: tau_ratqs
104	INTEGER :: iflag_cldcon, iflag_ratqs, iflag_coupl, iflag_clos, iflag_wake
105	INTEGER :: iflag_thermals, nsplit_thermals
106	INTEGER :: iflag_thermals_ed, iflag_thermals_optflux
107	REAL :: tau_thermals, solarlong0, seuil_inversion
108	INTEGER :: read_climoz ! read ozone climatology
109	! Allowed values are 0, 1 and 2
110	! 0: do not read an ozone climatology
111	! 1: read a single ozone climatology that will be used day and night
112	! 2: read two ozone climatologies, the average day and night
113	! climatology and the daylight climatology
114	!-------------------------------------------------------------------------------
115	REAL :: alp_offset
116	logical found
117
118	!--- Constants
119	pi = 4. * ATAN(1.)
120	rad = 6371229.
121	daysec = 86400.
122	omeg = 2.*pi/daysec
123	g = 9.8
124	kappa = 0.2857143
125	cpp = 1004.70885
126	preff = 101325.
127	pa = 50000.
128	jmp1 = jjm + 1
129
130	!--- CONSTRUCT A GRID
131	CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, &
132	callstats, &
133	solarlong0,seuil_inversion, &
134	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
135	iflag_cldcon, &
136	iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
137	ok_ade, ok_aie, aerosol_couple, &
138	flag_aerosol, new_aod, &
139	bl95_b0, bl95_b1, &
140	iflag_thermals,nsplit_thermals,tau_thermals, &
141	iflag_thermals_ed,iflag_thermals_optflux, &
142	iflag_coupl,iflag_clos,iflag_wake, read_climoz, &
143	alp_offset)
144
145	! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
146	co2_ppm0 = co2_ppm
147
148	dtvr = daysec/FLOAT(day_step)
149	WRITE(lunout,*)'dtvr',dtvr
150
151	CALL iniconst()
152	CALL inigeom()
153
154	!--- Initializations for tracers
155	CALL infotrac_init
156	ALLOCATE(q3d(iip1,jjp1,llm,nqtot))
157
158	CALL inifilr()
159	CALL phys_state_var_init(read_climoz)
160
161	rlat(1) = ASIN(1.0)
162	DO j=2,jjm; rlat((j-2)iim+2:(j-1)iim+1)=rlatu(j); END DO
163	rlat(klon) = - ASIN(1.0)
164	rlat(:)=rlat(:)*(180.0/pi)
165
166	rlon(1) = 0.0
167	DO j=2,jjm; rlon((j-2)iim+2:(j-1)iim+1)=rlonv(1:iim); END DO
168	rlon(klon) = 0.0
169	rlon(:)=rlon(:)*(180.0/pi)
170
171	! For a coupled simulation, the ocean mask from ocean model is used to compute
172	! the weights an to insure ocean fractions are the same for atmosphere and ocean
173	! Otherwise, mask is created using Relief file.
174
175	WRITE(lunout,*)'Essai de lecture masque ocean'
176	iret = NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)
177	IF(iret/=NF90_NOERR) THEN
178	WRITE(lunout,*)'ATTENTION!! pas de fichier o2a.nc trouve'
179	WRITE(lunout,*)'Run force'
180	x='masque'
181	masque(:,:)=0.0
182	CALL startget_phys2d(x, iip1, jjp1, rlonv, rlatu, masque, 0.0, jjm, &
183	& rlonu, rlatv, ib)
184	WRITE(lunout,*)'MASQUE construit : Masque'
185	WRITE(lunout,'(97I1)') nINT(masque)
186	CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
187	WHERE( zmasq(:)<EPSFRA) zmasq(:)=0.
188	WHERE(1.-zmasq(:)<EPSFRA) zmasq(:)=1.
189	ELSE
190	WRITE(lunout,*)'ATTENTION!! fichier o2a.nc trouve'
191	WRITE(lunout,*)'Run couple'
192	couple=.true.
193	iret=NF90_CLOSE(nid_o2a)
194	CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
195	IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
196	WRITE(lunout,*)'Dimensions non compatibles pour masque ocean'
197	WRITE(lunout,*)'iim = ',iim,' iml_omask = ',iml_omask
198	WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
199	CALL abort_gcm('etat0_netcdf','Dimensions non compatibles pour masque oc&
200	&ean',1)
201	END IF
202	ALLOCATE( ocemask(iml_omask,jml_omask), ocetmp(iml_omask,jml_omask))
203	ALLOCATE( lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask))
204	ALLOCATE(dlon_omask(iml_omask), dlat_omask(jml_omask))
205	CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, lon_omask,&
206	& lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
207	CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, ttm_tmp, &
208	& 1, 1, ocetmp)
209	CALL flinclo(fid)
210	dlon_omask(1:iml_omask) = lon_omask(1:iml_omask,1)
211	dlat_omask(1:jml_omask) = lat_omask(1,1:jml_omask)
212	ocemask = ocetmp
213	IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
214	DO j=1,jml_omask
215	ocemask(:,j) = ocetmp(:,jml_omask-j+1)
216	END DO
217	END IF
218	!
219	! Ocean mask to physical grid
220	!*******************************************************************************
221	WRITE(lunout,*)'ocemask '
222	WRITE(fmt,"(i4,'i1)')")iml_omask ; fmt='('//ADJUSTL(fmt)
223	WRITE(lunout,fmt)int(ocemask)
224	ocemask_fi(1)=ocemask(1,1)
225	DO j=2,jjm; ocemask_fi((j-2)iim+2:(j-1)iim+1)=ocemask(1:iim,j); END DO
226	ocemask_fi(klon)=ocemask(1,jjp1)
227	zmasq=1.-ocemask_fi
228	END IF
229
230	CALL gr_fi_dyn(1,klon,iip1,jjp1,zmasq,masque)
231
232	! The startget calls need to be replaced by a call to restget to get the
233	! values in the restart file
234	x = 'relief'; orog(:,:) = 0.0
235	CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu, orog, 0.0,jjm,rlonu,rlatv,ib,&
236	& masque)
237
238	x = 'rugosite'; rugo(:,:) = 0.0
239	CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu, rugo, 0.0,jjm, rlonu,rlatv,ib)
240	! WRITE(lunout,'(49I1)') INT(orog(:,:)*10)
241	! WRITE(lunout,'(49I1)') INT(rugo(:,:)*10)
242
243	! Sub-surfaces initialization
244	!*******************************************************************************
245	pctsrf=0.
246	x = 'psol'; psol(:,:) = 0.0
247	CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu,psol,0.0,jjm,rlonu,rlatv,ib)
248	! WRITE(lunout,*) 'PSOL :', psol(10,20)
249	! WRITE(lunout,*) ap(:), bp(:)
250
251	! Mid-levels pressure computation
252	!*******************************************************************************
253	CALL pression(ip1jmp1, ap, bp, psol, p3d)
254	CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
255	pls(:,:,:)=preff(pk(:,:,:)/cpp)*(1./kappa)
256	! WRITE(lunout,*) 'P3D :', p3d(10,20,:)
257	! WRITE(lunout,*) 'PK:', pk(10,20,:)
258	! WRITE(lunout,*) 'PLS :', pls(10,20,:)
259
260	x = 'surfgeo'; phis(:,:) = 0.0
261	CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu,phis, 0.0,jjm, rlonu,rlatv,ib)
262
263	x = 'u'; uvent(:,:,:) = 0.0
264	CALL startget_dyn(x,rlonu,rlatu,pls,y,uvent,0.0, &
265	& rlonv,rlatv,ib)
266
267	x = 'v'; vvent(:,:,:) = 0.0
268	CALL startget_dyn(x, rlonv,rlatv,pls(:, :jjm, :),y(:, :jjm, :),vvent,0.0, &
269	& rlonu,rlatu(:jjm),ib)
270
271	x = 't'; t3d(:,:,:) = 0.0
272	CALL startget_dyn(x,rlonv,rlatu,pls,y,t3d,0.0, &
273	& rlonu,rlatv,ib)
274
275	x = 'tpot'; tpot(:,:,:) = 0.0
276	CALL startget_dyn(x,rlonv,rlatu,pls,pk,tpot,0.0, &
277	& rlonu,rlatv,ib)
278
279	WRITE(lunout,*) 'T3D min,max:',minval(t3d(:,:,:)),maxval(t3d(:,:,:))
280	WRITE(lunout,*) 'PLS min,max:',minval(pls(:,:,:)),maxval(pls(:,:,:))
281
282	! Humidity at saturation computation
283	!*******************************************************************************
284	WRITE(lunout,*) 'avant q_sat'
285	CALL q_sat(llmjjp1iip1, t3d, pls, qsat)
286	WRITE(lunout,*) 'apres q_sat'
287	WRITE(lunout,*) 'QSAT min,max:',minval(qsat(:,:,:)),maxval(qsat(:,:,:))
288	! WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
289
290	x = 'q'; qd (:,:,:) = 0.0
291	CALL startget_dyn(x,rlonv,rlatu,pls,qsat,qd,0.0, rlonu,rlatv,ib)
292	q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:)
293
294	! Parameterization of ozone chemistry:
295	! Look for ozone tracer:
296	i = 1
297	DO
298	found = tname(i)=="O3" .OR. tname(i)=="o3"
299	if (found .or. i == nqtot) exit
300	i = i + 1
301	end do
302	if (found) then
303	call regr_lat_time_coefoz
304	call press_coefoz
305	call regr_pr_o3(p3d, q3d(:, :, :, i))
306	! Convert from mole fraction to mass fraction:
307	q3d(:, :, :, i) = q3d(:, :, :, i) * 48. / 29.
308	end if
309
310	!--- OZONE CLIMATOLOGY
311	IF(read_climoz>=1) CALL regr_lat_time_climoz(read_climoz)
312
313	x = 'tsol'; tsol(:) = 0.0
314	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,tsol,0.0,jjm,rlonu,rlatv,ib)
315
316	x = 'qsol'; qsol(:) = 0.0
317	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,qsol,0.0,jjm,rlonu,rlatv,ib)
318
319	x = 'snow'; sn(:) = 0.0
320	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,sn,0.0,jjm,rlonu,rlatv,ib)
321
322	x = 'rads'; radsol(:) = 0.0
323	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,jjm,rlonu,rlatv,ib)
324
325	x = 'rugmer'; rugmer(:) = 0.0
326	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,jjm,rlonu,rlatv,ib)
327
328	x = 'zmea'; zmea(:) = 0.0
329	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,jjm,rlonu,rlatv,ib)
330
331	x = 'zstd'; zstd(:) = 0.0
332	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,jjm,rlonu,rlatv,ib)
333
334	x = 'zsig'; zsig(:) = 0.0
335	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,jjm,rlonu,rlatv,ib)
336
337	x = 'zgam'; zgam(:) = 0.0
338	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,jjm,rlonu,rlatv,ib)
339
340	x = 'zthe'; zthe(:) = 0.0
341	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,jjm,rlonu,rlatv,ib)
342
343	x = 'zpic'; zpic(:) = 0.0
344	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,jjm,rlonu,rlatv,ib)
345
346	x = 'zval'; zval(:) = 0.0
347	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,jjm,rlonu,rlatv,ib)
348
349	! WRITE(lunout,'(48I3)') 'TSOL :', INT(tsol(2:klon)-273)
350
351	! Soil ice file reading for soil fraction and soil ice fraction
352	!*******************************************************************************
353	CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
354	ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
355	ALLOCATE(dlat_lic(jml_lic), dlon_lic(iml_lic))
356	ALLOCATE( fraclic(iml_lic,jml_lic))
357	CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp, &
358	& lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
359	CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
360	CALL flinclo(fid)
361
362	! Interpolation on model T-grid
363	!*******************************************************************************
364	WRITE(lunout,*)'dimensions de landice iml_lic, jml_lic : ',iml_lic,jml_lic
365	! conversion if coordinates are in degrees
366	IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180.
367	IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
368	dlon_lic(:)=lon_lic(:,1)
369	dlat_lic(:)=lat_lic(1,:)
370	CALL grille_m( iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic, iim,jjp1, &
371	& rlonv, rlatu, flic_tmp(1:iim,:) )
372	flic_tmp(iip1,:)=flic_tmp(1,:)
373
374	!--- To the physical grid
375	CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic))
376
377	!--- Adequation with soil/sea mask
378	WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0.
379	WHERE(zmasq(:)<EPSFRA) pctsrf(:,is_lic)=0.
380	pctsrf(:,is_ter)=zmasq(:)
381	DO ji=1,klon
382	IF(zmasq(ji)>EPSFRA) THEN
383	IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
384	pctsrf(ji,is_lic)=zmasq(ji)
385	pctsrf(ji,is_ter)=0.
386	ELSE
387	pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
388	IF(pctsrf(ji,is_ter)<EPSFRA) THEN
389	pctsrf(ji,is_ter)=0.
390	pctsrf(ji,is_lic)=zmasq(ji)
391	END IF
392	END IF
393	END IF
394	END DO
395
396	! sub-surface ocean and sea ice (sea ice set to zero for start)
397	!*******************************************************************************
398	pctsrf(:,is_oce)=(1.-zmasq(:))
399	WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
400	IF(couple) pctsrf(:,is_oce)=ocemask_fi(:)
401	isst=0
402	WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
403
404	! It is checked that the sub-surfaces sum is equal to 1
405	!*******************************************************************************
406	ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
407	IF(ji/=0) WRITE(lunout,*) 'pb repartition sous maille pour ',ji,' points'
408	CALL gr_fi_ecrit(1, klon, iim, jjp1, zmasq, zx_tmp_2d)
409	! WRITE(fmt,"(i3,')')")iim; fmt='(i'//ADJUSTL(fmt)
410	! WRITE(lunout,*)'zmasq = '
411	! WRITE(lunout,TRIM(fmt))NINT(zx_tmp_2d)
412	CALL gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
413	WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
414	WRITE(lunout,*) 'MASQUE construit : Masque'
415	WRITE(lunout,TRIM(fmt))NINT(masque(:,:))
416
417	! Intermediate computation
418	!*******************************************************************************
419	CALL massdair(p3d,masse)
420	WRITE(lunout,*)' ALPHAX ',alphax
421	DO l=1,llm
422	xppn(:)=aire(1:iim,1 )*masse(1:iim,1 ,l)
423	xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
424	xpn=SUM(xppn)/apoln
425	xps=SUM(xpps)/apols
426	masse(:,1 ,l)=xpn
427	masse(:,jjp1,l)=xps
428	END DO
429	q3d(iip1,:,:,:)=q3d(1,:,:,:)
430	phis(iip1,:) = phis(1,:)
431
432	IF(.NOT.letat0) RETURN
433
434	! Writing
435	!*******************************************************************************
436	CALL inidissip(lstardis,nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,tetatemp)
437	WRITE(lunout,*)'sortie inidissip'
438	itau=0
439	itau_dyn=0
440	itau_phy=0
441	iday=dayref+itau/day_step
442	time=FLOAT(itau-(iday-dayref)*day_step)/day_step
443	IF(time>1.) THEN
444	time=time-1
445	iday=iday+1
446	END IF
447	day_ref=dayref
448	annee_ref=anneeref
449
450	CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
451	WRITE(lunout,*)'sortie geopot'
452
453	CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, &
454	phi, w, pbaru, pbarv, time+iday-dayref)
455	WRITE(lunout,*)'sortie caldyn0'
456	CALL dynredem0( "start.nc", dayref, phis)
457	WRITE(lunout,*)'sortie dynredem0'
458	CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
459	WRITE(lunout,*)'sortie dynredem1'
460
461	! Physical initial state writting
462	!*******************************************************************************
463	WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
464	phystep = dtvr * FLOAT(iphysiq)
465	radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) )
466	WRITE(lunout,*)'phystep =', phystep, radpas
467
468	! Init: tsol, qsol, sn, evap, tsoil, rain_fall, snow_fall, solsw, sollw, frugs
469	!*******************************************************************************
470	DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
471	DO i=1,nbsrf; snsrf(:,i) = sn; END DO
472	falb1(:,is_ter) = 0.08; falb1(:,is_lic) = 0.6
473	falb1(:,is_oce) = 0.5; falb1(:,is_sic) = 0.6
474	falb2 = falb1
475	evap(:,:) = 0.
476	DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
477	DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
478	rain_fall = 0.; snow_fall = 0.
479	solsw = 165.; sollw = -53.
480	t_ancien = 273.15
481	q_ancien = 0.
482	agesno = 0.
483	frugs(:,is_oce) = rugmer(:)
484	frugs(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
485	frugs(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
486	frugs(:,is_sic) = 0.001
487	fder = 0.0
488	clwcon = 0.0
489	rnebcon = 0.0
490	ratqs = 0.0
491	run_off_lic_0 = 0.0
492	rugoro = 0.0
493
494	! Before phyredem calling, surface modules and values to be saved in startphy.nc
495	! are initialized
496	!*******************************************************************************
497	dummy = 1.0
498	pbl_tke(:,:,:) = 1.e-8
499	zmax0(:) = 40.
500	f0(:) = 1.e-5
501	ema_work1(:,:) = 0.
502	ema_work2(:,:) = 0.
503	wake_deltat(:,:) = 0.
504	wake_deltaq(:,:) = 0.
505	wake_s(:) = 0.
506	wake_cstar(:) = 0.
507	wake_fip(:) = 0.
508	wake_pe = 0.
509	fm_therm = 0.
510	entr_therm = 0.
511	detr_therm = 0.
512
513	CALL fonte_neige_init(run_off_lic_0)
514	CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil )
515	CALL phyredem( "startphy.nc" )
516
517	! WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
518	! WRITE(lunout,*)'entree histclo'
519	CALL histclo()
520
521	#endif
522	!#endif of #ifdef CPP_EARTH
523	RETURN
524
525	END SUBROUTINE etat0_netcdf
526	!
527	!-------------------------------------------------------------------------------

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

Download in other formats:

Original Format