Context Navigation

etat0_netcdf.F @ 5444

Last change on this file since 5444 was 546, checked in by lmdzadmin, 20 years ago
Suite des adaptations au portage sur SGI et VPP pour Prism (Caubel & Demory) LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.1 KB

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

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

Download in other formats:

Original Format