Context Navigation

etat0_netcdf.F90 @ 1320

Last change on this file since 1320 was 1319, checked in by Laurent Fairhead, 15 years ago

Modifications to the start and limit creation routines to account for different

calendars

Modification to phyetat0 to force the mask read in the start files to match the

surface fractions read in the limit file

Force readaerosol.F90 to read in aerosols file with 12 timesteps

Modifications aux routines de créations des fichiers start et limit pour prendre

en compte différents calendriers

Modification à phyetat0 pour forcer le masque lu dans le fichier start à être

compatible avec les fractions de surface lu dans le fichier limit

Forcer readaerosol à ne lire que des fichiers à 12 pas de temps

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

File size: 19.0 KB

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

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

Download in other formats:

Original Format