Context Navigation

etat0_netcdf.F @ 1176

Last change on this file since 1176 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 23.7 KB

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

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

Download in other formats:

Original Format