Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

etat0_netcdf.F @ 1368

Last change on this file since 1368 was 701, checked in by (none), 19 years ago
This commit was manufactured by cvs2svn to create branch 'pre_V3'.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.1 KB

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

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

Download in other formats:

Original Format