Context Navigation

etat0dyn_netcdf.F90 @ 2600

Last change on this file since 2600 was 2600, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90 EM
File size: 21.7 KB

Line
1	MODULE etat0dyn
2	!
3	!*******************************************************************************
4	! Purpose: Create dynamical initial state using atmospheric fields from a
5	! database of atmospheric to initialize the model.
6	!-------------------------------------------------------------------------------
7	! Comments:
8	!
9	! * This module is designed to work for Earth (and with ioipsl)
10	!
11	! * etat0dyn_netcdf routine can access to NetCDF data through the following
12	! routine (to be called after restget):
13	! CALL startget_dyn3d(varname, lon_in, lat_in, pls, workvar,&
14	! champ, lon_in2, lat_in2)
15	!
16	! * Variables should have the following names in the NetCDF files:
17	! 'U' : East ward wind (in "ECDYN.nc")
18	! 'V' : Northward wind (in "ECDYN.nc")
19	! 'TEMP' : Temperature (in "ECDYN.nc")
20	! 'R' : Relative humidity (in "ECDYN.nc")
21	! 'RELIEF' : High resolution orography (in "Relief.nc")
22	!
23	! * The land mask and corresponding weights can be:
24	! 1) already known (in particular if etat0dyn has been called before) ;
25	! in this case, ANY(masque(:,:)/=-99999.) = .TRUE.
26	! 2) computed using the ocean mask from the ocean model (to ensure ocean
27	! fractions are the same for atmosphere and ocean) for coupled runs.
28	! File name: "o2a.nc" ; variable name: "OceMask"
29	! 3) computed from topography file "Relief.nc" for forced runs.
30	!
31	! * There is a big mess with the longitude size. Should it be iml or iml+1 ?
32	! I have chosen to use the iml+1 as an argument to this routine and we declare
33	! internaly smaller fields when needed. This needs to be cleared once and for
34	! all in LMDZ. A convention is required.
35	!-------------------------------------------------------------------------------
36	USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, histclo
37	USE assert_eq_m, ONLY: assert_eq
38	USE comconst_mod, ONLY: pi, cpp, kappa
39	USE comvert_mod, ONLY: ap, bp, preff, pressure_exner
40	IMPLICIT NONE
41
42	PRIVATE
43	PUBLIC :: etat0dyn_netcdf
44
45	include "iniprint.h"
46	include "dimensions.h"
47	include "paramet.h"
48	include "comgeom2.h"
49	include "temps.h"
50	include "comdissnew.h"
51	REAL, SAVE :: deg2rad
52	INTEGER, SAVE :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
53	REAL, ALLOCATABLE, SAVE :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)
54	CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc'
55
56	CONTAINS
57
58	!-------------------------------------------------------------------------------
59	!
60	SUBROUTINE etat0dyn_netcdf(masque, phis)
61	!
62	!-------------------------------------------------------------------------------
63	! Purpose: Create dynamical initial states.
64	!-------------------------------------------------------------------------------
65	! Notes: 1) This routine is designed to work for Earth
66	! 2) If masque(:,:)/=-99999., masque and phis are already known.
67	! Otherwise: compute it.
68	!-------------------------------------------------------------------------------
69	USE control_mod
70	USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
71	USE regr_pr_o3_m, ONLY: regr_pr_o3
72	USE press_coefoz_m, ONLY: press_coefoz
73	USE exner_hyb_m, ONLY: exner_hyb
74	USE exner_milieu_m, ONLY: exner_milieu
75	USE infotrac, ONLY: nqtot, tname
76	USE filtreg_mod
77	IMPLICIT NONE
78	!-------------------------------------------------------------------------------
79	! Arguments:
80	REAL, INTENT(INOUT) :: masque(iip1,jjp1) !--- Land-ocean mask
81	REAL, INTENT(INOUT) :: phis (iip1,jjp1) !--- Ground geopotential
82	!-------------------------------------------------------------------------------
83	! Local variables:
84	CHARACTER(LEN=256) :: modname, fmt
85	INTEGER :: i, j, l, ji, itau, iday
86	REAL :: xpn, xps, time, phystep
87	REAL, DIMENSION(iip1,jjp1) :: psol
88	REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
89	REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd
90	REAL, DIMENSION(iip1,jjp1,llm) :: pk, pls, y, masse
91	REAL, DIMENSION(iip1,jjm ,llm) :: vvent
92	REAL, DIMENSION(ip1jm ,llm) :: pbarv
93	REAL, DIMENSION(ip1jmp1 ,llm) :: pbaru, phi, w
94	REAL, DIMENSION(ip1jmp1) :: pks
95	REAL, DIMENSION(iim) :: xppn, xpps
96	REAL, ALLOCATABLE :: q3d(:,:,:,:)
97	!-------------------------------------------------------------------------------
98	modname='etat0dyn_netcdf'
99
100	deg2rad = pi/180.0
101
102	! Compute psol AND tsol, knowing phis.
103	!*******************************************************************************
104	CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol)
105
106	! Mid-levels pressure computation
107	!*******************************************************************************
108	CALL pression(ip1jmp1, ap, bp, psol, p3d) !--- Update p3d
109	IF(pressure_exner) THEN !--- Update pk, pks
110	CALL exner_hyb (ip1jmp1,psol,p3d,pks,pk)
111	ELSE
112	CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk)
113	END IF
114	pls(:,:,:)=preff(pk(:,:,:)/cpp)*(1./kappa) !--- Update pls
115
116	! Update uvent, vvent, t3d and tpot
117	!*******************************************************************************
118	uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
119	CALL startget_dyn3d('u' ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv)
120	CALL startget_dyn3d('v' ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent, &
121	& rlonu,rlatu(:jjm))
122	CALL startget_dyn3d('t' ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv)
123	tpot(:,:,:)=t3d(:,:,:)
124	CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv)
125
126	WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
127	WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:))
128
129	! Humidity at saturation computation
130	!*******************************************************************************
131	WRITE(lunout,*) 'avant q_sat'
132	CALL q_sat(llmjjp1iip1, t3d, pls, qsat)
133	WRITE(lunout,*) 'apres q_sat'
134	WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:))
135	! WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
136	qd (:,:,:) = 0.0
137	CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv)
138	ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
139	CALL flinclo(fid_dyn)
140
141	! Parameterization of ozone chemistry:
142	!*******************************************************************************
143	! Look for ozone tracer:
144	DO i=1,nqtot; IF(ANY(["O3","o3"]==tname(i))) EXIT; END DO
145	IF(i/=nqtot+1) THEN
146	CALL regr_lat_time_coefoz
147	CALL press_coefoz
148	CALL regr_pr_o3(p3d, q3d(:,:,:,i))
149	q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29. !--- Mole->mass fraction
150	END IF
151	q3d(iip1,:,:,:)=q3d(1,:,:,:)
152
153	! Writing
154	!*******************************************************************************
155	CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, &
156	tetatemp, vert_prof_dissip)
157	WRITE(lunout,*)'sortie inidissip'
158	itau=0
159	itau_dyn=0
160	itau_phy=0
161	iday=dayref+itau/day_step
162	time=FLOAT(itau-(iday-dayref)*day_step)/day_step
163	IF(time>1.) THEN
164	time=time-1
165	iday=iday+1
166	END IF
167	day_ref=dayref
168	annee_ref=anneeref
169	CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
170	WRITE(lunout,*)'sortie geopot'
171	CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, &
172	phi, w, pbaru, pbarv, time+iday-dayref)
173	WRITE(lunout,*)'sortie caldyn0'
174	#ifdef CPP_PARA
175	CALL dynredem0_loc( "start.nc", dayref, phis)
176	#else
177	CALL dynredem0( "start.nc", dayref, phis)
178	#endif
179	WRITE(lunout,*)'sortie dynredem0'
180	#ifdef CPP_PARA
181	CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
182	#else
183	CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
184	#endif
185	WRITE(lunout,*)'sortie dynredem1'
186	CALL histclo()
187
188	END SUBROUTINE etat0dyn_netcdf
189	!
190	!-------------------------------------------------------------------------------
191
192
193	!-------------------------------------------------------------------------------
194	!
195	SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar,&
196	champ, lon_in2, lat_in2)
197	!-------------------------------------------------------------------------------
198	IMPLICIT NONE
199	!===============================================================================
200	! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R
201	! (3D fields) of file dynfname.
202	!-------------------------------------------------------------------------------
203	! Note: An input auxilliary field "workvar" has to be specified in two cases:
204	! * for "q": the saturated humidity.
205	! * for "tpot": the Exner function.
206	!===============================================================================
207	! Arguments:
208	CHARACTER(LEN=*), INTENT(IN) :: var
209	REAL, INTENT(IN) :: lon_in(:) ! dim (iml)
210	REAL, INTENT(IN) :: lat_in(:) ! dim (jml)
211	REAL, INTENT(IN) :: pls (:, :, :) ! dim (iml, jml, lml)
212	REAL, INTENT(IN) :: workvar(:, :, :) ! dim (iml, jml, lml)
213	REAL, INTENT(INOUT) :: champ (:, :, :) ! dim (iml, jml, lml)
214	REAL, INTENT(IN) :: lon_in2(:) ! dim (iml)
215	REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2)
216	!-------------------------------------------------------------------------------
217	! Local variables:
218	CHARACTER(LEN=10) :: vname
219	CHARACTER(LEN=256) :: msg, modname="startget_dyn3d"
220	INTEGER :: iml, jml, jml2, lml, il
221	REAL :: xppn, xpps
222	!-------------------------------------------------------------------------------
223	iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1), &
224	& SIZE(lon_in2)], TRIM(modname)//" iml")
225	jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2), &
226	& TRIM(modname)//" jml")
227	lml=assert_eq( SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3), &
228	& TRIM(modname)//" lml")
229	jml2=SIZE(lat_in2)
230
231	!--- CHECK IF THE FIELD IS KNOWN
232	SELECT CASE(var)
233	CASE('u'); vname='U'
234	CASE('v'); vname='V'
235	CASE('t'); vname='TEMP'
236	CASE('q'); vname='R'; msg='humidity as the saturated humidity'
237	CASE('tpot'); msg='potential temperature as the Exner function'
238	CASE DEFAULT; msg='No rule to extract variable '//TRIM(var)
239	CALL abort_gcm(modname,TRIM(msg)//' from any data set',1)
240	END SELECT
241
242	!--- CHECK IF SOMETHING IS MISSING
243	IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
244	msg='Could not compute '//TRIM(msg)//' is missing or constant.'
245	CALL abort_gcm(modname,TRIM(msg),1)
246	END IF
247
248	!--- INTERPOLATE 3D FIELD IF NEEDED
249	IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2, &
250	lat_in2,pls,champ)
251
252	!--- COMPUTE THE REQUIRED FILED
253	SELECT CASE(var)
254	CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
255	champ(iml,:,:)=champ(1,:,:) !--- Eastward wind
256
257	CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
258	champ(iml,:,:)=champ(1,:,:) !--- Northward wind
259
260	CASE('tpot','q')
261	IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature
262	ELSE; champ=champ.01workvar !--- Relative humidity
263	WHERE(champ<0.) champ=1.0E-10
264	END IF
265	DO il=1,lml
266	xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln
267	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
268	champ(:,1 ,il) = xppn
269	champ(:,jml,il) = xpps
270	END DO
271	END SELECT
272
273	END SUBROUTINE startget_dyn3d
274	!
275	!-------------------------------------------------------------------------------
276
277
278	!-------------------------------------------------------------------------------
279	!
280	SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol)
281	!
282	!-------------------------------------------------------------------------------
283	IMPLICIT NONE
284	!===============================================================================
285	! Purpose: Compute psol, knowing phis.
286	!===============================================================================
287	! Arguments:
288	REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml)
289	REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2)
290	REAL, INTENT(IN) :: zs (:,:) ! dim (iml,jml)
291	REAL, INTENT(OUT) :: psol(:,:) ! dim (iml,jml)
292	!-------------------------------------------------------------------------------
293	! Local variables:
294	CHARACTER(LEN=256) :: modname='start_init_dyn'
295	REAL :: date, dt
296	INTEGER :: iml, jml, jml2, itau(1)
297	REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:,:)
298	REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:)
299	REAL, ALLOCATABLE :: z(:,:), ps(:,:), ts(:,:)
300	!-------------------------------------------------------------------------------
301	iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2), &
302	& TRIM(modname)//" iml")
303	jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml")
304	jml2=SIZE(lat_in2)
305
306	WRITE(lunout,*) 'Opening the surface analysis'
307	CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
308	WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
309
310	ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn))
311	ALLOCATE(levdyn_ini(llm_dyn))
312	CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, &
313	lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn)
314
315	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
316	ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
317	lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
318	lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
319
320	ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
321	CALL get_var_dyn('Z',z) !--- SURFACE GEOPOTENTIAL
322	CALL get_var_dyn('SP',ps) !--- SURFACE PRESSURE
323	CALL get_var_dyn('ST',ts) !--- SURFACE TEMPERATURE
324	! CALL flinclo(fid_dyn)
325	DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
326
327	!--- PSOL IS COMPUTED IN PASCALS
328	psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0 &
329	& /ts(:iml-1,:))
330	psol(iml,:)=psol(1,:)
331	DEALLOCATE(z,ps,ts)
332	psol(:,1 )=SUM(aire(1:iml-1,1 )*psol(1:iml-1,1 ))/apoln !--- NORTH POLE
333	psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols !--- SOUTH POLE
334
335	CONTAINS
336
337	!-------------------------------------------------------------------------------
338	!
339	SUBROUTINE get_var_dyn(title,field)
340	!
341	!-------------------------------------------------------------------------------
342	USE conf_dat_m, ONLY: conf_dat2d
343	IMPLICIT NONE
344	!-------------------------------------------------------------------------------
345	! Arguments:
346	CHARACTER(LEN=*), INTENT(IN) :: title
347	REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
348	!-------------------------------------------------------------------------------
349	! Local variables:
350	CHARACTER(LEN=256) :: msg
351	INTEGER :: tllm
352	!-------------------------------------------------------------------------------
353	SELECT CASE(title)
354	CASE('Z'); tllm=0; msg='geopotential'
355	CASE('SP'); tllm=0; msg='surface pressure'
356	CASE('ST'); tllm=llm_dyn; msg='temperature'
357	END SELECT
358	IF(.NOT.ALLOCATED(field)) THEN
359	ALLOCATE(field(iml,jml))
360	CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
361	CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
362	CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana, &
363	lon_in, lat_in, lon_in2, lat_in2, field)
364	ELSE IF(SIZE(field)/=SIZE(z)) THEN
365	msg='The '//TRIM(msg)//' field we have does not have the right size'
366	CALL abort_gcm(TRIM(modname),msg,1)
367	END IF
368
369	END SUBROUTINE get_var_dyn
370	!
371	!-------------------------------------------------------------------------------
372
373	END SUBROUTINE start_init_dyn
374	!
375	!-------------------------------------------------------------------------------
376
377
378	!-------------------------------------------------------------------------------
379	!
380	SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d)
381	!
382	!-------------------------------------------------------------------------------
383	USE conf_dat_m, ONLY: conf_dat3d
384	USE pchsp_95_m, ONLY: pchsp_95
385	USE pchfe_95_m, ONLY: pchfe_95
386	IMPLICIT NONE
387	!-------------------------------------------------------------------------------
388	! Arguments:
389	CHARACTER(LEN=*), INTENT(IN) :: var
390	REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml)
391	REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2)
392	REAL, INTENT(IN) :: pls_in(:,:,:) ! dim (iml,jml,lml)
393	REAL, INTENT(OUT) :: var3d (:,:,:) ! dim (iml,jml,lml)
394	!-------------------------------------------------------------------------------
395	! Local variables:
396	CHARACTER(LEN=256) :: modname='start_inter_3d'
397	LOGICAL :: skip
398	REAL :: chmin, chmax
399	INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr
400	INTEGER :: n_extrap ! Extrapolated points number
401	REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:)
402	REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:)
403	REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:)
404	!-------------------------------------------------------------------------------
405	iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml")
406	jml=assert_eq(SIZE(lat_in), SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml")
407	lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2)
408
409	WRITE(lunout, *)'Going into flinget to extract the 3D field.'
410	IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
411	CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
412
413	!--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS
414	ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
415	lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
416	lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
417
418	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
419	ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn))
420	CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, &
421	lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.)
422	DEALLOCATE(lon_ini, lat_ini)
423
424	!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro
425	ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
426	DO il = 1,llm_dyn
427	CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il), &
428	lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
429	END DO
430	DEALLOCATE(lon_rad, lat_rad)
431
432	!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND
433	ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
434	ax = lev_dyn(llm_dyn:1:-1)
435	skip = .FALSE.
436	n_extrap = 0
437	DO ij=1, jml
438	DO ii=1, iml-1
439	ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
440	yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
441	CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
442	var3d(ii, ij, lml:1:-1), ierr)
443	IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1)
444	n_extrap = n_extrap + ierr
445	END DO
446	END DO
447	IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap
448	var3d(iml, :, :) = var3d(1, :, :)
449
450	DO il=1, lml
451	CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
452	WRITE(lunout, *)' '//TRIM(var)//' min max l ', il, chmin, chmax
453	END DO
454
455	END SUBROUTINE start_inter_3d
456	!
457	!-------------------------------------------------------------------------------
458
459
460	!-------------------------------------------------------------------------------
461	!
462	SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
463	!
464	!-------------------------------------------------------------------------------
465	USE inter_barxy_m, ONLY: inter_barxy
466	IMPLICIT NONE
467	!-------------------------------------------------------------------------------
468	! Arguments:
469	CHARACTER(LEN=*), INTENT(IN) :: nam
470	LOGICAL, INTENT(IN) :: ibeg
471	REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj)
472	REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj)
473	REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1)
474	REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2)
475	REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1)
476	!-------------------------------------------------------------------------------
477	! Local variables:
478	CHARACTER(LEN=256) :: modname="interp_startvar"
479	INTEGER :: ii, jj, i1, j1, j2
480	REAL, ALLOCATABLE :: vtmp(:,:)
481	!-------------------------------------------------------------------------------
482	ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")
483	jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")
484	i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
485	j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1")
486	j2=SIZE(lat2)
487	ALLOCATE(vtmp(i1-1,j1))
488	IF(ibeg.AND.prt_level>1) THEN
489	WRITE(lunout,*)"---------------------------------------------------------"
490	WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
491	WRITE(lunout,*)"---------------------------------------------------------"
492	END IF
493	CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
494	CALL gr_int_dyn(vtmp, varo, i1-1, j1)
495
496	END SUBROUTINE interp_startvar
497	!
498	!-------------------------------------------------------------------------------
499
500	END MODULE etat0dyn
501	!
502	!*******************************************************************************

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

Download in other formats:

Original Format