Context Navigation

etat0_netcdf.F @ 1212

Last change on this file since 1212 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.6 KB

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

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

Download in other formats:

Original Format