Context Navigation

etat0dyn_netcdf.F90 @ 5281

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

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

Download in other formats:

Original Format