Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

startvar.F @ 5416

Last change on this file since 5416 was 253, checked in by (none), 23 years ago
This commit was manufactured by cvs2svn to create branch 'rel-1-0-patch'.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 32.6 KB

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

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

Download in other formats:

Original Format