Context Navigation

startvar.F @ 273

Last change on this file since 273 was 259, checked in by lmdz, 23 years ago
Nouveaux programmes pour la creation des etats initiaux et des conditions aux limites. Levan LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 38.9 KB

Rev	Line
[259]	1	C
	2	C $Header$
	3	C
[2]	4	MODULE startvar
	5	!
	6	!
	7	! There are three ways to access data from the database of atmospheric data which
	8	! can be used to initialize the model. This depends on the type of field which needs
	9	! to be extracted. In any case the call should come after a restget and should be of the type :
	10	! CALL startget(...)
	11	!
	12	! We will details the possible arguments to startget here :
	13	!
	14	! - A 2D variable on the dynamical grid :
[259]	15	! CALL startget(varname, iml, jml, lon_in, lat_in, champ, val_ex, jml2, lon_in2, lat_in2, interbar )
[2]	16	!
	17	! - A 1D variable on the physical grid :
[259]	18	! CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
[2]	19	!
	20	!
	21	! - A 3D variable on the dynamical grid :
[259]	22	! CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp, jml2, lon_in2, lat_in2, interbar )
[2]	23	!
	24	!
	25	! There is special constraint on the atmospheric data base except that the
	26	! the data needs to be in netCDF and the variables should have the the following
	27	! names in the file :
	28	!
	29	! 'RELIEF' : High resolution orography
	30	! 'ST' : Surface temperature
	31	! 'CDSW' : Soil moisture
	32	! 'Z' : Surface geopotential
	33	! 'SP' : Surface pressure
	34	! 'U' : East ward wind
	35	! 'V' : Northward wind
	36	! 'TEMP' : Temperature
	37	! 'R' : Relative humidity
	38	!
	39	USE ioipsl
	40	!
	41	!
	42	IMPLICIT NONE
	43	!
	44	!
	45	PRIVATE
	46	PUBLIC startget
	47	!
	48	!
	49	INTERFACE startget
	50	MODULE PROCEDURE startget_phys2d, startget_phys1d, startget_dyn
	51	END INTERFACE
	52	!
	53	INTEGER, SAVE :: fid_phys, fid_dyn
	54	INTEGER, SAVE :: iml_phys, iml_rel, iml_dyn
	55	INTEGER, SAVE :: jml_phys, jml_rel, jml_dyn
	56	INTEGER, SAVE :: llm_dyn, ttm_dyn
	57	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: lon_phys, lon_rug,
	58	. lon_alb, lon_rel, lon_dyn
	59	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: lat_phys, lat_rug,
	60	. lat_alb, lat_rel, lat_dyn
[259]	61	REAL, ALLOCATABLE, SAVE, DIMENSION (:) :: levdyn_ini
[2]	62	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: relief, zstd, zsig,
	63	. zgam, zthe, zpic, zval
	64	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rugo, masque, phis
	65	!
	66	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: tsol, qsol, psol_dyn
	67	REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: var_ana3d
	68	!
	69	CONTAINS
	70	!
	71	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	72	!
	73	SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in,
[259]	74	. champ, val_exp, jml2, lon_in2, lat_in2 , interbar )
[2]	75	!
	76	! There is a big mess with the size in logitude, should it be iml or iml+1.
	77	! I have chosen to use the iml+1 as an argument to this routine and we declare
	78	! internaly smaler fields when needed. This needs to be cleared once and for all in LMDZ.
	79	! A convention is required.
	80	!
	81	!
	82	CHARACTER(), INTENT(in) :: varname
[259]	83	INTEGER, INTENT(in) :: iml, jml ,jml2
[2]	84	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
[259]	85	REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
[2]	86	REAL, INTENT(inout) :: champ(iml,jml)
	87	REAL, INTENT(in) :: val_exp
[259]	88	LOGICAL interbar
[2]	89	!
	90	! This routine only works if the variable does not exist or is constant
	91	!
	92	IF ( MINVAL(champ(:,:)).EQ.MAXVAL(champ(:,:)) .AND.
	93	.MINVAL(champ(:,:)).EQ.val_exp ) THEN
	94	!
	95	SELECTCASE(varname)
	96	!
	97	CASE ('relief')
	98	!
	99	! If we do not have the orography we need to get it
	100	!
	101	IF ( .NOT.ALLOCATED(relief)) THEN
	102	!
[259]	103	CALL start_init_orog( iml, jml, lon_in, lat_in,
	104	. jml2,lon_in2,lat_in2, interbar )
[2]	105	!
	106	ENDIF
	107	!
	108	IF (SIZE(relief) .NE. SIZE(champ)) THEN
	109	!
	110	WRITE(,) 'STARTVAR module has been',
	111	.' initialized to the wrong size'
	112	STOP
	113	!
	114	ENDIF
	115	!
	116	champ(:,:) = relief(:,:)
	117	!
	118	CASE ('rugosite')
	119	!
	120	! If we do not have the orography we need to get it
	121	!
	122	IF ( .NOT.ALLOCATED(rugo)) THEN
	123	!
[259]	124	CALL start_init_orog( iml, jml, lon_in, lat_in,
	125	. jml2,lon_in2,lat_in2 , interbar )
[2]	126	!
	127	ENDIF
	128	!
	129	IF (SIZE(rugo) .NE. SIZE(champ)) THEN
	130	!
	131	WRITE(,)
	132	. 'STARTVAR module has been initialized to the wrong size'
	133	STOP
	134	!
	135	ENDIF
	136	!
	137	champ(:,:) = rugo(:,:)
	138	!
	139	CASE ('masque')
	140	!
	141	! If we do not have the orography we need to get it
	142	!
	143	IF ( .NOT.ALLOCATED(masque)) THEN
	144	!
[259]	145	CALL start_init_orog( iml, jml, lon_in, lat_in,
	146	. jml2,lon_in2,lat_in2 , interbar )
[2]	147	!
	148	ENDIF
	149	!
	150	IF (SIZE(masque) .NE. SIZE(champ)) THEN
	151	!
	152	WRITE(,)
	153	. 'STARTVAR module has been initialized to the wrong size'
	154	STOP
	155	!
	156	ENDIF
	157	!
	158	champ(:,:) = masque(:,:)
	159	!
	160	CASE ('surfgeo')
	161	!
	162	! If we do not have the orography we need to get it
	163	!
	164	IF ( .NOT.ALLOCATED(phis)) THEN
	165	!
[259]	166	CALL start_init_orog( iml, jml, lon_in, lat_in,
	167	. jml2,lon_in2, lat_in2 , interbar )
[2]	168	!
	169	ENDIF
	170	!
	171	IF (SIZE(phis) .NE. SIZE(champ)) THEN
	172	!
	173	WRITE(,)
	174	. 'STARTVAR module has been initialized to the wrong size'
	175	STOP
	176	!
	177	ENDIF
	178	!
	179	champ(:,:) = phis(:,:)
	180	!
	181	CASE ('psol')
	182	!
	183	! If we do not have the orography we need to get it
	184	!
	185	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
	186	!
[259]	187	CALL start_init_dyn( iml, jml, lon_in, lat_in,
	188	. jml2,lon_in2, lat_in2 , interbar )
[2]	189	!
	190	ENDIF
	191	!
	192	IF (SIZE(psol_dyn) .NE. SIZE(champ)) THEN
	193	!
	194	WRITE(,)
	195	. 'STARTVAR module has been initialized to the wrong size'
	196	STOP
	197	!
	198	ENDIF
	199	!
	200	champ(:,:) = psol_dyn(:,:)
	201	!
	202	CASE DEFAULT
	203	!
	204	WRITE(,) 'startget_phys2d'
	205	WRITE(,) 'No rule is present to extract variable',
	206	. varname(:LEN_TRIM(varname)),' from any data set'
	207	STOP
	208	!
	209	END SELECT
	210	!
	211	ELSE
	212	!
	213	! There are a few fields we might need if we need to interpolate 3D filed. Thus if they come through here we
	214	! will catch them
	215	!
	216	SELECTCASE(varname)
	217	!
	218	CASE ('surfgeo')
	219	!
	220	IF ( .NOT.ALLOCATED(phis)) THEN
	221	ALLOCATE(phis(iml,jml))
	222	ENDIF
	223	!
	224	IF (SIZE(phis) .NE. SIZE(champ)) THEN
	225	!
	226	WRITE(,)
	227	. 'STARTVAR module has been initialized to the wrong size'
	228	STOP
	229	!
	230	ENDIF
	231	!
	232	phis(:,:) = champ(:,:)
	233	!
	234	END SELECT
	235	!
	236	ENDIF
	237	!
	238	END SUBROUTINE startget_phys2d
	239	!
	240	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	241	!
[259]	242	SUBROUTINE start_init_orog ( iml,jml,lon_in, lat_in,jml2,lon_in2 ,
	243	, lat_in2 , interbar )
[2]	244	!
[259]	245	INTEGER, INTENT(in) :: iml, jml, jml2
[2]	246	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
[259]	247	REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
	248	LOGICAL interbar
[2]	249	!
	250	! LOCAL
	251	!
[259]	252	LOGICAL interbar2
	253	REAL :: lev(1), date, dt,chmin,chmax
[2]	254	INTEGER :: itau(1), fid
	255	INTEGER :: llm_tmp, ttm_tmp
	256	INTEGER :: i, j
	257	INTEGER :: iret
[259]	258	CHARACTER*25 title
[2]	259	REAL, ALLOCATABLE :: relief_hi(:,:)
	260	REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
[259]	261	REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:)
[2]	262	REAL, ALLOCATABLE :: tmp_var(:,:)
	263	INTEGER, ALLOCATABLE :: tmp_int(:,:)
	264	!
	265	CHARACTER*120 :: orogfname
	266	LOGICAL :: check=.TRUE.
	267	!
	268	!
	269	orogfname = 'Relief.nc'
	270	!
	271	IF ( check ) WRITE(,) 'Reading the high resolution orography'
	272	!
	273	CALL flininfo(orogfname,iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
	274	!
	275	ALLOCATE (lat_rel(iml_rel,jml_rel), stat=iret)
	276	ALLOCATE (lon_rel(iml_rel,jml_rel), stat=iret)
	277	ALLOCATE (relief_hi(iml_rel,jml_rel), stat=iret)
	278	!
	279	CALL flinopen(orogfname, .FALSE., iml_rel, jml_rel,
	280	.llm_tmp, lon_rel, lat_rel, lev, ttm_tmp,
	281	. itau, date, dt, fid)
	282	!
	283	CALL flinget(fid, 'RELIEF', iml_rel, jml_rel, llm_tmp,
	284	. ttm_tmp, 1, 1, relief_hi)
	285	!
	286	CALL flinclo(fid)
	287	!
	288	! In case we have a file which is in degrees we do the transformation
	289	!
	290	ALLOCATE(lon_rad(iml_rel))
[259]	291	ALLOCATE(lon_ini(iml_rel))
	292
[2]	293	IF ( MAXVAL(lon_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
[259]	294	lon_ini(:) = lon_rel(:,1) * 2.0 * ASIN(1.0) / 180.0
[2]	295	ELSE
[259]	296	lon_ini(:) = lon_rel(:,1)
[2]	297	ENDIF
[259]	298
[2]	299	ALLOCATE(lat_rad(jml_rel))
[259]	300	ALLOCATE(lat_ini(jml_rel))
	301
[2]	302	IF ( MAXVAL(lat_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
[259]	303	lat_ini(:) = lat_rel(1,:) * 2.0 * ASIN(1.0) / 180.0
[2]	304	ELSE
[259]	305	lat_ini(:) = lat_rel(1,:)
[2]	306	ENDIF
	307	!
	308	!
[259]	309
	310	title='RELIEF'
	311
	312	interbar2 = .FALSE.
	313	CALL conf_dat2d(title,iml_rel, jml_rel, lon_ini, lat_ini,
	314	. lon_rad, lat_rad, relief_hi , interbar2 )
	315
[2]	316	IF ( check ) WRITE(,) 'Computes all the parameters needed',
	317	.' for the gravity wave drag code'
	318	!
	319	! Allocate the data we need to put in the interpolated fields
	320	!
	321	! RELIEF: orographie moyenne
	322	ALLOCATE(relief(iml,jml))
	323	! zphi : orographie moyenne
	324	ALLOCATE(phis(iml,jml))
	325	! zstd: deviation standard de l'orographie sous-maille
	326	ALLOCATE(zstd(iml,jml))
	327	! zsig: pente de l'orographie sous-maille
	328	ALLOCATE(zsig(iml,jml))
	329	! zgam: anisotropy de l'orographie sous maille
	330	ALLOCATE(zgam(iml,jml))
	331	! zthe: orientation de l'axe oriente dans la direction
	332	! de plus grande pente de l'orographie sous maille
	333	ALLOCATE(zthe(iml,jml))
	334	! zpic: hauteur pics de la SSO
	335	ALLOCATE(zpic(iml,jml))
	336	! zval: hauteur vallees de la SSO
	337	ALLOCATE(zval(iml,jml))
	338	! masque : Masque terre ocean
	339	ALLOCATE(tmp_int(iml,jml))
	340	ALLOCATE(masque(iml,jml))
	341	!
	342	CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi,
	343	. iml-1, jml, lon_in, lat_in,
	344	. phis, relief, zstd, zsig, zgam, zthe, zpic, zval, tmp_int)
[24]	345	phis = phis * 9.81
[2]	346	!
	347	masque(:,:) = FLOAT(tmp_int(:,:))
	348	!
	349	! Compute surface roughness
	350	!
	351	IF ( check ) WRITE(,)
	352	.'Compute surface roughness induced by the orography'
	353	!
	354	ALLOCATE(rugo(iml,jml))
	355	ALLOCATE(tmp_var(iml-1,jml))
	356	!
	357	CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi,
	358	. iml-1, jml, lon_in, lat_in, tmp_var)
	359	!
	360	DO j = 1, jml
	361	DO i = 1, iml-1
	362	rugo(i,j) = tmp_var(i,j)
	363	ENDDO
	364	rugo(iml,j) = tmp_var(1,j)
	365	ENDDO
[259]	366	c
	367	cc * rugo n'est pas utilise pour l'instant ****
[2]	368	!
	369	! Build land-sea mask
	370	!
	371	!
	372	RETURN
	373	!
	374	END SUBROUTINE start_init_orog
	375	!
	376	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	377	!
	378	SUBROUTINE startget_phys1d(varname, iml, jml, lon_in,
[259]	379	.lat_in, nbindex, champ, val_exp ,jml2, lon_in2, lat_in2,interbar)
[2]	380	!
	381	CHARACTER(), INTENT(in) :: varname
[259]	382	INTEGER, INTENT(in) :: iml, jml, nbindex, jml2
[2]	383	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
[259]	384	REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
[2]	385	REAL, INTENT(inout) :: champ(nbindex)
	386	REAL, INTENT(in) :: val_exp
[259]	387	LOGICAL interbar
[2]	388	!
	389	!
	390	! This routine only works if the variable does not exist or is constant
	391	!
	392	IF ( MINVAL(champ(:)).EQ.MAXVAL(champ(:)) .AND.
	393	.MINVAL(champ(:)).EQ.val_exp ) THEN
	394	SELECTCASE(varname)
	395	CASE ('tsol')
	396	IF ( .NOT.ALLOCATED(tsol)) THEN
[259]	397	CALL start_init_phys( iml, jml, lon_in, lat_in,
	398	. jml2, lon_in2, lat_in2, interbar )
[2]	399	ENDIF
	400	IF ( SIZE(tsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
	401	WRITE(,)
	402	. 'STARTVAR module has been initialized to the wrong size'
	403	STOP
	404	ENDIF
	405	CALL gr_dyn_fi(1, iml, jml, nbindex, tsol, champ)
	406	CASE ('qsol')
	407	IF ( .NOT.ALLOCATED(qsol)) THEN
[259]	408	CALL start_init_phys( iml, jml, lon_in, lat_in,
	409	. jml2, lon_in2,lat_in2 , interbar )
[2]	410	ENDIF
	411	IF ( SIZE(qsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
	412	WRITE(,)
	413	. 'STARTVAR module has been initialized to the wrong size'
	414	STOP
	415	ENDIF
	416	CALL gr_dyn_fi(1, iml, jml, nbindex, qsol, champ)
	417	CASE ('psol')
	418	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
[259]	419	CALL start_init_dyn( iml, jml, lon_in, lat_in,
	420	. jml2, lon_in2,lat_in2 , interbar )
[2]	421	ENDIF
	422	IF (SIZE(psol_dyn) .NE. SIZE(lon_in)*SIZE(lat_in)) THEN
	423	WRITE(,)
	424	. 'STARTVAR module has been initialized to the wrong size'
	425	STOP
	426	ENDIF
	427	CALL gr_dyn_fi(1, iml, jml, nbindex, psol_dyn, champ)
	428	CASE ('zmea')
	429	IF ( .NOT.ALLOCATED(relief)) THEN
[259]	430	CALL start_init_orog( iml, jml, lon_in, lat_in,
	431	. jml2, lon_in2,lat_in2 , interbar )
[2]	432	ENDIF
	433	IF ( SIZE(relief) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
	434	WRITE(,)
	435	. 'STARTVAR module has been initialized to the wrong size'
	436	STOP
	437	ENDIF
	438	CALL gr_dyn_fi(1, iml, jml, nbindex, relief, champ)
	439	CASE ('zstd')
	440	IF ( .NOT.ALLOCATED(zstd)) THEN
[259]	441	CALL start_init_orog( iml, jml, lon_in, lat_in,
	442	. jml2, lon_in2,lat_in2 , interbar )
[2]	443	ENDIF
	444	IF ( SIZE(zstd) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
	445	WRITE(,)
	446	. 'STARTVAR module has been initialized to the wrong size'
	447	STOP
	448	ENDIF
	449	CALL gr_dyn_fi(1, iml, jml, nbindex,zstd, champ)
	450	CASE ('zsig')
	451	IF ( .NOT.ALLOCATED(zsig)) THEN
[259]	452	CALL start_init_orog( iml, jml, lon_in, lat_in,
	453	. jml2, lon_in2,lat_in2 , interbar )
[2]	454	ENDIF
	455	IF ( SIZE(zsig) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
	456	WRITE(,)
	457	. 'STARTVAR module has been initialized to the wrong size'
	458	STOP
	459	ENDIF
	460	CALL gr_dyn_fi(1, iml, jml, nbindex,zsig, champ)
	461	CASE ('zgam')
	462	IF ( .NOT.ALLOCATED(zgam)) THEN
[259]	463	CALL start_init_orog( iml, jml, lon_in, lat_in,
	464	. jml2, lon_in2,lat_in2 , interbar )
[2]	465	ENDIF
	466	IF ( SIZE(zgam) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
	467	WRITE(,)
	468	. 'STARTVAR module has been initialized to the wrong size'
	469	STOP
	470	ENDIF
	471	CALL gr_dyn_fi(1, iml, jml, nbindex,zgam, champ)
	472	CASE ('zthe')
	473	IF ( .NOT.ALLOCATED(zthe)) THEN
[259]	474	CALL start_init_orog( iml, jml, lon_in, lat_in,
	475	. jml2, lon_in2,lat_in2 , interbar )
[2]	476	ENDIF
	477	IF ( SIZE(zthe) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
	478	WRITE(,)
	479	. 'STARTVAR module has been initialized to the wrong size'
	480	STOP
	481	ENDIF
	482	CALL gr_dyn_fi(1, iml, jml, nbindex,zthe, champ)
	483	CASE ('zpic')
	484	IF ( .NOT.ALLOCATED(zpic)) THEN
[259]	485	CALL start_init_orog( iml, jml, lon_in, lat_in,
	486	. jml2, lon_in2,lat_in2 , interbar )
[2]	487	ENDIF
	488	IF ( SIZE(zpic) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
	489	WRITE(,)
	490	. 'STARTVAR module has been initialized to the wrong size'
	491	STOP
	492	ENDIF
	493	CALL gr_dyn_fi(1, iml, jml, nbindex,zpic, champ)
	494	CASE ('zval')
	495	IF ( .NOT.ALLOCATED(zval)) THEN
[259]	496	CALL start_init_orog( iml, jml, lon_in, lat_in,
	497	. jml2, lon_in2,lat_in2 , interbar )
[2]	498	ENDIF
	499	IF ( SIZE(zval) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
	500	WRITE(,)
	501	. 'STARTVAR module has been initialized to the wrong size'
	502	STOP
	503	ENDIF
	504	CALL gr_dyn_fi(1, iml, jml, nbindex,zval, champ)
[259]	505	CASE ('rads')
[2]	506	champ(:) = 0.0
	507	CASE ('snow')
	508	champ(:) = 0.0
	509	CASE ('deltat')
	510	champ(:) = 0.0
	511	CASE ('rugmer')
	512	champ(:) = 0.001
	513	CASE ('agsno')
	514	champ(:) = 50.0
	515	CASE DEFAULT
	516	WRITE(,) 'startget_phys1d'
	517	WRITE(,) 'No rule is present to extract variable ',
	518	. varname(:LEN_TRIM(varname)),' from any data set'
	519	STOP
	520	END SELECT
	521	ELSE
	522	!
	523	! If we see tsol we catch it as we may need it for a 3D interpolation
	524	!
	525	SELECTCASE(varname)
	526	CASE ('tsol')
	527	IF ( .NOT.ALLOCATED(tsol)) THEN
	528	ALLOCATE(tsol(SIZE(lon_in),SIZE(lat_in) ))
	529	ENDIF
	530	CALL gr_fi_dyn(1, iml, jml, nbindex, champ, tsol)
	531	END SELECT
	532	ENDIF
	533	END SUBROUTINE startget_phys1d
	534	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	535	!
[259]	536	SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in, jml2,
	537	. lon_in2, lat_in2 , interbar )
[2]	538	!
[259]	539	INTEGER, INTENT(in) :: iml, jml ,jml2
[2]	540	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
[259]	541	REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
	542	LOGICAL interbar
[2]	543	!
	544	! LOCAL
	545	!
	546	REAL :: lev(1), date, dt
	547	INTEGER :: itau(1)
	548	INTEGER :: llm_tmp, ttm_tmp
	549	INTEGER :: i, j
	550	!
[259]	551	CHARACTER*25 title
[2]	552	CHARACTER*120 :: physfname
	553	LOGICAL :: check=.TRUE.
	554	!
	555	REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
[259]	556	REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:)
[2]	557	REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:)
	558	!
	559	physfname = 'ECPHY.nc'
	560	!
	561	IF ( check ) WRITE(,) 'Opening the surface analysis'
	562	!
	563	CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp,
	564	. ttm_tmp, fid_phys)
	565	!
	566	!
	567	ALLOCATE (lat_phys(iml_phys,jml_phys))
	568	ALLOCATE (lon_phys(iml_phys,jml_phys))
	569	!
	570	CALL flinopen(physfname, .FALSE., iml_phys, jml_phys,
	571	. llm_tmp, lon_phys, lat_phys, lev, ttm_tmp,
	572	. itau, date, dt, fid_phys)
	573	!
	574	! Allocate the space we will need to get the data out of this file
	575	!
	576	ALLOCATE(var_ana(iml_phys, jml_phys))
	577	!
	578	! In case we have a file which is in degrees we do the transformation
	579	!
	580	ALLOCATE(lon_rad(iml_phys))
[259]	581	ALLOCATE(lon_ini(iml_phys))
	582
[2]	583	IF ( MAXVAL(lon_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
[259]	584	lon_ini(:) = lon_phys(:,1) * 2.0 * ASIN(1.0) / 180.0
[2]	585	ELSE
[259]	586	lon_ini(:) = lon_phys(:,1)
[2]	587	ENDIF
[259]	588
[2]	589	ALLOCATE(lat_rad(jml_phys))
[259]	590	ALLOCATE(lat_ini(jml_phys))
	591
[2]	592	IF ( MAXVAL(lat_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
[259]	593	lat_ini(:) = lat_phys(1,:) * 2.0 * ASIN(1.0) / 180.0
[2]	594	ELSE
[259]	595	lat_ini(:) = lat_phys(1,:)
[2]	596	ENDIF
[259]	597
	598
[2]	599	!
	600	! We get the two standard varibales
	601	! Surface temperature
	602	!
	603	ALLOCATE(tsol(iml,jml))
	604	ALLOCATE(tmp_var(iml-1,jml))
	605	!
	606	!
[259]	607
[2]	608	CALL flinget(fid_phys, 'ST', iml_phys, jml_phys,
	609	.llm_tmp, ttm_tmp, 1, 1, var_ana)
[259]	610
	611	title='ST'
	612	CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini,
	613	. lon_rad, lat_rad, var_ana , interbar )
	614
	615	IF ( interbar ) THEN
	616	WRITE(6,*) '-------------------------------------------------',
	617	,'--------------'
	618	WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
	619	, ' pour ST $$$ '
	620	WRITE(6,*) '-------------------------------------------------',
	621	,'--------------'
	622	CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad ,
	623	, var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var )
	624	ELSE
	625	CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
	626	. var_ana, iml-1, jml, lon_in, lat_in, tmp_var )
	627	ENDIF
	628
[2]	629	CALL gr_int_dyn(tmp_var, tsol, iml-1, jml)
	630	!
	631	! Soil moisture
	632	!
	633	ALLOCATE(qsol(iml,jml))
	634	CALL flinget(fid_phys, 'CDSW', iml_phys, jml_phys,
	635	. llm_tmp, ttm_tmp, 1, 1, var_ana)
[259]	636
	637	title='CDSW'
	638	CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini,
	639	. lon_rad, lat_rad, var_ana, interbar )
	640
	641	IF ( interbar ) THEN
	642	WRITE(6,*) '-------------------------------------------------',
	643	,'--------------'
	644	WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
	645	, ' pour CDSW $$$ '
	646	WRITE(6,*) '-------------------------------------------------',
	647	,'--------------'
	648	CALL inter_barxy ( iml_phys,jml_phys -1,lon_rad,lat_rad ,
	649	, var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var )
	650	ELSE
	651	CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
	652	. var_ana, iml-1, jml, lon_in, lat_in, tmp_var )
	653	ENDIF
	654	c
	655	CALL gr_int_dyn(tmp_var, qsol, iml-1, jml)
[2]	656	!
[259]	657	CALL flinclo(fid_phys)
[2]	658	!
	659	END SUBROUTINE start_init_phys
	660	!
	661	!
	662	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	663	!
	664	!
	665	SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in,
[259]	666	. lml, pls, workvar, champ, val_exp,jml2, lon_in2, lat_in2 ,
	667	, interbar )
[2]	668	!
	669	! ARGUMENTS
	670	!
	671	CHARACTER(), INTENT(in) :: varname
[259]	672	INTEGER, INTENT(in) :: iml, jml, lml, jml2
[2]	673	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
[259]	674	REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
[2]	675	REAL, INTENT(in) :: pls(iml, jml, lml)
	676	REAL, INTENT(in) :: workvar(iml, jml, lml)
	677	REAL, INTENT(inout) :: champ(iml, jml, lml)
	678	REAL, INTENT(in) :: val_exp
[259]	679	LOGICAL interbar
[2]	680	!
	681	! LOCAL
	682	!
	683	INTEGER :: il, ij, ii
	684	REAL :: xppn, xpps
	685	!
	686	#include "dimensions.h"
	687	#include "paramet.h"
	688	#include "comgeom2.h"
	689	#include "comconst.h"
	690	!
	691	! This routine only works if the variable does not exist or is constant
	692	!
	693	IF ( MINVAL(champ(:,:,:)).EQ.MAXVAL(champ(:,:,:)) .AND.
	694	. MINVAL(champ(:,:,:)).EQ.val_exp ) THEN
	695	!
	696	SELECTCASE(varname)
	697	CASE ('u')
	698	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
[259]	699	CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 ,
	700	. lon_in2,lat_in2 , interbar )
[2]	701	ENDIF
	702	CALL start_inter_3d('U', iml, jml, lml, lon_in,
[259]	703	. lat_in, jml2, lon_in2, lat_in2, pls, champ,interbar )
[2]	704	DO il=1,lml
	705	DO ij=1,jml
	706	DO ii=1,iml-1
	707	champ(ii,ij,il) = champ(ii,ij,il) * cu(ii,ij)
	708	ENDDO
	709	champ(iml,ij, il) = champ(1,ij, il)
	710	ENDDO
	711	ENDDO
	712	CASE ('v')
	713	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
[259]	714	CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2,
	715	. lon_in2, lat_in2 , interbar )
[2]	716	ENDIF
[259]	717	CALL start_inter_3d('V', iml, jml, lml, lon_in,
	718	. lat_in, jml2, lon_in2, lat_in2, pls, champ, interbar )
[2]	719	DO il=1,lml
	720	DO ij=1,jml
	721	DO ii=1,iml-1
	722	champ(ii,ij,il) = champ(ii,ij,il) * cv(ii,ij)
	723	ENDDO
	724	champ(iml,ij, il) = champ(1,ij, il)
	725	ENDDO
	726	ENDDO
	727	CASE ('t')
	728	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
[259]	729	CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 ,
	730	. lon_in2, lat_in2 ,interbar )
[2]	731	ENDIF
	732	CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
[259]	733	. lat_in, jml2, lon_in2, lat_in2, pls, champ, interbar )
[2]	734
	735	CASE ('tpot')
	736	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
[259]	737	CALL start_init_dyn( iml, jml, lon_in, lat_in , jml2 ,
	738	. lon_in2, lat_in2 , interbar )
[2]	739	ENDIF
	740	CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
[259]	741	. lat_in, jml2, lon_in2, lat_in2, pls, champ, interbar )
[2]	742	IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
	743	. THEN
	744	DO il=1,lml
	745	DO ij=1,jml
	746	DO ii=1,iml-1
	747	champ(ii,ij,il) = champ(ii,ij,il) * cpp
	748	. / workvar(ii,ij,il)
	749	ENDDO
	750	champ(iml,ij,il) = champ(1,ij,il)
	751	ENDDO
	752	ENDDO
	753	DO il=1,lml
	754	xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
[19]	755	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
[2]	756	champ(:,1,il) = xppn
	757	champ(:,jml,il) = xpps
	758	ENDDO
	759	ELSE
	760	WRITE(,)'Could not compute potential temperature as the'
	761	WRITE(,)'Exner function is missing or constant.'
	762	STOP
	763	ENDIF
	764	CASE ('q')
	765	IF ( .NOT.ALLOCATED(psol_dyn)) THEN
[259]	766	CALL start_init_dyn( iml, jml, lon_in, lat_in, jml2 ,
	767	. lon_in2, lat_in2 , interbar )
[2]	768	ENDIF
	769	CALL start_inter_3d('R', iml, jml, lml, lon_in, lat_in,
[259]	770	. jml2, lon_in2, lat_in2, pls, champ, interbar )
[2]	771	IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
	772	. THEN
	773	DO il=1,lml
	774	DO ij=1,jml
	775	DO ii=1,iml-1
	776	champ(ii,ij,il) = 0.01 * champ(ii,ij,il) *
	777	. workvar(ii,ij,il)
	778	ENDDO
	779	champ(iml,ij,il) = champ(1,ij,il)
	780	ENDDO
	781	ENDDO
	782	WHERE ( champ .LT. 0.) champ = 1.0E-10
	783	DO il=1,lml
	784	xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
[19]	785	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
[2]	786	champ(:,1,il) = xppn
	787	champ(:,jml,il) = xpps
	788	ENDDO
	789	ELSE
	790	WRITE(,)'Could not compute specific humidity as the'
	791	WRITE(,)'saturated humidity is missing or constant.'
	792	STOP
	793	ENDIF
	794	CASE DEFAULT
	795	WRITE(,) 'startget_dyn'
	796	WRITE(,) 'No rule is present to extract variable ',
	797	. varname(:LEN_TRIM(varname)),' from any data set'
	798	STOP
	799	END SELECT
	800	ENDIF
	801	END SUBROUTINE startget_dyn
	802	!
	803	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	804	!
[259]	805	SUBROUTINE start_init_dyn( iml, jml, lon_in, lat_in,jml2,lon_in2 ,
	806	, lat_in2 , interbar )
[2]	807	!
[259]	808	INTEGER, INTENT(in) :: iml, jml, jml2
[2]	809	REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
[259]	810	REAL, INTENT(in) :: lon_in2(iml), lat_in2(jml2)
	811	LOGICAL interbar
[2]	812	!
	813	! LOCAL
	814	!
	815	REAL :: lev(1), date, dt
	816	INTEGER :: itau(1)
	817	INTEGER :: i, j
	818	integer :: iret
	819	!
	820	CHARACTER*120 :: physfname
	821	LOGICAL :: check=.TRUE.
	822	!
	823	REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
[259]	824	REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:)
[2]	825	REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:), z(:,:)
	826	REAL, ALLOCATABLE :: xppn(:), xpps(:)
	827	LOGICAL :: allo
	828	!
	829	!
	830	#include "dimensions.h"
	831	#include "paramet.h"
	832	#include "comgeom2.h"
[259]	833
	834	CHARACTER*25 title
	835
[2]	836	!
	837	physfname = 'ECDYN.nc'
	838	!
	839	IF ( check ) WRITE(,) 'Opening the surface analysis'
	840	!
	841	CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn,
	842	. ttm_dyn, fid_dyn)
	843	IF ( check ) WRITE(,) 'Values read: ', iml_dyn, jml_dyn,
	844	. llm_dyn, ttm_dyn
	845	!
	846	ALLOCATE (lat_dyn(iml_dyn,jml_dyn), stat=iret)
	847	ALLOCATE (lon_dyn(iml_dyn,jml_dyn), stat=iret)
[259]	848	ALLOCATE (levdyn_ini(llm_dyn), stat=iret)
[2]	849	!
	850	CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,
[259]	851	. lon_dyn, lat_dyn, levdyn_ini, ttm_dyn,
[2]	852	. itau, date, dt, fid_dyn)
	853	!
	854
	855	allo = allocated (var_ana)
	856	if (allo) then
	857	DEALLOCATE(var_ana, stat=iret)
	858	endif
	859	ALLOCATE(var_ana(iml_dyn, jml_dyn), stat=iret)
	860
	861	allo = allocated (lon_rad)
	862	if (allo) then
	863	DEALLOCATE(lon_rad, stat=iret)
	864	endif
[259]	865
	866	ALLOCATE(lon_rad(iml_dyn), stat=iret)
	867	ALLOCATE(lon_ini(iml_dyn))
[2]	868
	869	IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
[259]	870	lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
[2]	871	ELSE
[259]	872	lon_ini(:) = lon_dyn(:,1)
[2]	873	ENDIF
[259]	874
[2]	875	ALLOCATE(lat_rad(jml_dyn))
[259]	876	ALLOCATE(lat_ini(jml_dyn))
	877
[2]	878	IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
[259]	879	lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
[2]	880	ELSE
[259]	881	lat_ini(:) = lat_dyn(1,:)
[2]	882	ENDIF
	883	!
[259]	884
	885
[2]	886	ALLOCATE(z(iml, jml))
	887	ALLOCATE(tmp_var(iml-1,jml))
	888	!
	889	CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn,
	890	. 1, 1, var_ana)
[259]	891	c
	892	title='Z'
	893	CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini,
	894	. lon_rad, lat_rad, var_ana, interbar )
	895	c
	896	IF ( interbar ) THEN
	897	WRITE(6,*) '-------------------------------------------------',
	898	,'--------------'
	899	WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
	900	, ' pour Z $$$ '
	901	WRITE(6,*) '-------------------------------------------------',
	902	,'--------------'
	903	CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad ,
	904	, var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var)
	905	ELSE
	906	CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
[2]	907	. iml-1, jml, lon_in, lat_in, tmp_var)
[259]	908	ENDIF
	909
[2]	910	CALL gr_int_dyn(tmp_var, z, iml-1, jml)
	911	!
	912	ALLOCATE(psol_dyn(iml, jml))
	913	!
	914	CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn,
	915	. 1, 1, var_ana)
[259]	916
	917	title='SP'
	918	CALL conf_dat2d( title,iml_dyn, jml_dyn,lon_ini, lat_ini,
	919	. lon_rad, lat_rad, var_ana, interbar )
	920
	921	IF ( interbar ) THEN
	922	WRITE(6,*) '-------------------------------------------------',
	923	,'--------------'
	924	WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
	925	, ' pour SP $$$ '
	926	WRITE(6,*) '-------------------------------------------------',
	927	,'--------------'
	928	CALL inter_barxy ( iml_dyn,jml_dyn -1,lon_rad,lat_rad ,
	929	, var_ana, iml-1, jml-1, lon_in2, lat_in2, jml, tmp_var)
	930	ELSE
	931	CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
	932	. iml-1, jml, lon_in, lat_in, tmp_var )
	933	ENDIF
	934
[2]	935	CALL gr_int_dyn(tmp_var, psol_dyn, iml-1, jml)
	936	!
	937	IF ( .NOT.ALLOCATED(tsol)) THEN
	938	! These variables may have been allocated by the need to
	939	! create a start field for them or by the varibale
	940	! coming out of the restart file. In case we dor have it we will initialize it.
	941	!
[259]	942	CALL start_init_phys( iml, jml, lon_in, lat_in,jml2,lon_in2,
	943	. lat_in2 , interbar )
[2]	944	ELSE
	945	IF ( SIZE(tsol) .NE. SIZE(psol_dyn) ) THEN
	946	WRITE(,) 'start_init_dyn :'
	947	WRITE(,) 'The temperature field we have does not ',
	948	. 'have the right size'
	949	STOP
	950	ENDIF
	951	ENDIF
	952	IF ( .NOT.ALLOCATED(phis)) THEN
	953	!
	954	! These variables may have been allocated by the need to create a start field for them or by the varibale
	955	! coming out of the restart file. In case we dor have it we will initialize it.
	956	!
[259]	957	CALL start_init_orog( iml, jml, lon_in, lat_in, jml2, lon_in2 ,
	958	. lat_in2 , interbar )
[2]	959	!
	960	ELSE
	961	!
	962	IF (SIZE(phis) .NE. SIZE(psol_dyn)) THEN
	963	!
	964	WRITE(,) 'start_init_dyn :'
	965	WRITE(,) 'The orography field we have does not ',
	966	. ' have the right size'
	967	STOP
	968	ENDIF
	969	!
	970	ENDIF
	971	!
	972	! PSOL is computed in Pascals
	973	!
	974	!
	975	DO j = 1, jml
	976	DO i = 1, iml-1
	977	psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))
	978	. /287.0/tsol(i,j))
	979	ENDDO
	980	psol_dyn(iml,j) = psol_dyn(1,j)
	981	ENDDO
	982	!
	983	!
	984	ALLOCATE(xppn(iml-1))
	985	ALLOCATE(xpps(iml-1))
	986	!
	987	DO i = 1, iml-1
	988	xppn(i) = aire( i,1) * psol_dyn( i,1)
	989	xpps(i) = aire( i,jml) * psol_dyn( i,jml)
	990	ENDDO
	991	!
	992	DO i = 1, iml
	993	psol_dyn(i,1 ) = SUM(xppn)/apoln
	994	psol_dyn(i,jml) = SUM(xpps)/apols
	995	ENDDO
	996	!
	997	RETURN
	998	!
	999	END SUBROUTINE start_init_dyn
	1000	!
	1001	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1002	!
	1003	SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in,
[259]	1004	. lat_in, jml2, lon_in2, lat_in2, pls_in, var3d, interbar )
[2]	1005	!
	1006	! This subroutine gets a variables from a 3D file and does the interpolations needed
	1007	!
	1008	!
	1009	! ARGUMENTS
	1010	!
	1011	CHARACTER() :: varname
[259]	1012	INTEGER :: iml, jml, lml, jml2
[2]	1013	REAL :: lon_in(iml), lat_in(jml), pls_in(iml, jml, lml)
[259]	1014	REAL :: lon_in2(iml) , lat_in2(jml2)
[2]	1015	REAL :: var3d(iml, jml, lml)
[259]	1016	LOGICAL interbar
	1017	real chmin,chmax
[2]	1018	!
	1019	! LOCAL
	1020	!
[259]	1021	CHARACTER*25 title
	1022	INTEGER :: ii, ij, il, jsort,i,j,l
[2]	1023	REAL :: bx, by
	1024	REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
[259]	1025	REAL, ALLOCATABLE :: lon_ini(:), lat_ini(:) , lev_dyn(:)
[2]	1026	REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:)
	1027	REAL, ALLOCATABLE :: ax(:), ay(:), yder(:)
[259]	1028	REAL, ALLOCATABLE :: varrr(:,:,:)
[2]	1029	INTEGER, ALLOCATABLE :: lind(:)
	1030	!
	1031	LOGICAL :: check = .TRUE.
	1032	!
	1033	IF ( .NOT. ALLOCATED(var_ana3d)) THEN
	1034	ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
	1035	ENDIF
[259]	1036	ALLOCATE(varrr(iml_dyn, jml_dyn, llm_dyn))
[2]	1037	!
	1038	!
	1039	IF ( check) WRITE(,) 'Going into flinget to extract the 3D ',
	1040	. ' field.', fid_dyn
	1041	IF ( check) WRITE(,) fid_dyn, varname, iml_dyn, jml_dyn,
	1042	. llm_dyn,ttm_dyn
	1043	!
	1044	CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn,
	1045	. ttm_dyn, 1, 1, var_ana3d)
	1046	!
	1047	IF ( check) WRITE(,) 'Allocating space for the interpolation',
	1048	. iml, jml, llm_dyn
	1049	!
	1050	ALLOCATE(lon_rad(iml_dyn))
[259]	1051	ALLOCATE(lon_ini(iml_dyn))
	1052
[2]	1053	IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
[259]	1054	lon_ini(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
[2]	1055	ELSE
[259]	1056	lon_ini(:) = lon_dyn(:,1)
[2]	1057	ENDIF
[259]	1058
[2]	1059	ALLOCATE(lat_rad(jml_dyn))
[259]	1060	ALLOCATE(lat_ini(jml_dyn))
	1061
	1062	ALLOCATE(lev_dyn(llm_dyn))
	1063
[2]	1064	IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
[259]	1065	lat_ini(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
[2]	1066	ELSE
[259]	1067	lat_ini(:) = lat_dyn(1,:)
[2]	1068	ENDIF
	1069	!
[259]	1070
	1071	CALL conf_dat3d ( varname,iml_dyn, jml_dyn, llm_dyn, lon_ini,
	1072	. lat_ini, levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d ,
	1073	, interbar )
	1074
[2]	1075	ALLOCATE(var_tmp2d(iml-1, jml))
	1076	ALLOCATE(var_tmp3d(iml, jml, llm_dyn))
	1077	ALLOCATE(ax(llm_dyn))
	1078	ALLOCATE(ay(llm_dyn))
	1079	ALLOCATE(yder(llm_dyn))
	1080	ALLOCATE(lind(llm_dyn))
	1081	!
[259]	1082
[2]	1083	DO il=1,llm_dyn
	1084	!
[259]	1085	IF( interbar ) THEN
	1086	IF( il.EQ.1 ) THEN
	1087	WRITE(6,*) '-------------------------------------------------',
	1088	,'--------------'
	1089	WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
	1090	, ' pour ', varname
	1091	WRITE(6,*) '-------------------------------------------------',
	1092	,'--------------'
	1093	ENDIF
	1094	CALL inter_barxy ( iml_dyn, jml_dyn -1,lon_rad, lat_rad,
	1095	, var_ana3d(:,:,il),iml-1, jml2, lon_in2, lat_in2,jml,var_tmp2d )
	1096	ELSE
	1097	CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad,
	1098	. var_ana3d(:,:,il), iml-1, jml, lon_in, lat_in, var_tmp2d )
	1099	ENDIF
[2]	1100	!
	1101	CALL gr_int_dyn(var_tmp2d, var_tmp3d(:,:,il), iml-1, jml)
	1102	!
	1103	ENDDO
	1104	!
	1105	DO il=1,llm_dyn
	1106	lind(il) = llm_dyn-il+1
	1107	ENDDO
	1108	!
[259]	1109	c
	1110	c ... Pour l'interpolation verticale ,on interpole du haut de l'atmosphere
	1111	c vers le sol ...
	1112	c
[2]	1113	DO ij=1,jml
	1114	DO ii=1,iml-1
	1115	!
[259]	1116	ax(:) = lev_dyn(lind(:))
[2]	1117	ay(:) = var_tmp3d(ii, ij, lind(:))
	1118	!
[259]	1119
[2]	1120	CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
	1121	!
	1122	DO il=1,lml
	1123	bx = pls_in(ii, ij, il)
	1124	CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
	1125	var3d(ii, ij, il) = by
	1126	ENDDO
	1127	!
	1128	ENDDO
	1129	var3d(iml, ij, :) = var3d(1, ij, :)
	1130	ENDDO
	1131
[259]	1132	do il=1,lml
	1133	call minmax(iml*jml,var3d(1,1,il),chmin,chmax)
	1134	SELECTCASE(varname)
	1135	CASE('U')
	1136	WRITE(,) ' U min max l ',il,chmin,chmax
	1137	CASE('V')
	1138	WRITE(,) ' V min max l ',il,chmin,chmax
	1139	CASE('TEMP')
	1140	WRITE(,) ' TEMP min max l ',il,chmin,chmax
	1141	CASE('R')
	1142	WRITE(,) ' R min max l ',il,chmin,chmax
	1143	END SELECT
	1144	enddo
	1145
[2]	1146	DEALLOCATE(lon_rad)
	1147	DEALLOCATE(lat_rad)
	1148	DEALLOCATE(var_tmp2d)
	1149	DEALLOCATE(var_tmp3d)
	1150	DEALLOCATE(ax)
	1151	DEALLOCATE(ay)
	1152	DEALLOCATE(yder)
	1153	DEALLOCATE(lind)
	1154
	1155	!
	1156	RETURN
	1157	!
	1158	END SUBROUTINE start_inter_3d
	1159	!
	1160	END MODULE startvar

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

Download in other formats:

Original Format