Context Navigation

etat0_netcdf.F @ 1322

Last change on this file since 1322 was 1322, checked in by Laurent Fairhead, 15 years ago

Improvements concerning wake parametrisation (from JYG, NR, IT, with more to come).
Alp_offset is read in form physiq.def file

Améliorations à la paramétrisation des poches froides (de JYG, NR, IT, d'autres
sont à venir)
Alp_offset est rajouté à la liste des paramètres lus dans physiq.def

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 24.7 KB

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

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

Download in other formats:

Original Format