Context Navigation

etat0_netcdf.F90 @ 2239

Last change on this file since 2239 was 2239, checked in by Ehouarn Millour, 10 years ago

Reorganizing physics/dynamics interface:

what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
"bibio" directory is now "misc" and should only contain autonomous utilities.
"cosp" is now a subdirectory of phylmd.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

File size: 19.8 KB

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

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

Download in other formats:

Original Format