Context Navigation

etat0dyn_netcdf.F90 @ 2408

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

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

Download in other formats:

Original Format